Loading packages
#===============================================================================
#BTC.LineZero.Header.1.1.0
#===============================================================================
#R Markdown environment setup and reporting utility.
#===============================================================================
#RLB.Dependencies:
# knitr, magrittr, pacman, rio, rmarkdown, rmdformats, tibble, yaml
#===============================================================================
#Input for document parameters, libraries, file paths, and options.
#=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=+=
knitr::opts_chunk$set(message=FALSE, warning = FALSE)
path_working <- "/Users/minsikkim/Dropbox (Partners HealthCare)/Project_SICAS2_microbiome/5_Scripts/MGK/Host_depletion_git/"
path_library <- "/Library/Frameworks/R.framework/Resources/library"
str_libraries <- c(
"readxl", "phyloseq", "tidyverse", "pacman", "yaml"
)
path_working <- "/Users/minsikkim/Dropbox (Partners HealthCare)/Project_SICAS2_microbiome/5_Scripts/MGK/Host_depletion_git"
path_library <- "/Library/Frameworks/R.framework/Resources/library"
str_libraries <- c("readxl", "phyloseq", "tidyverse", "pacman", "yaml", "ggplot2", "vegan", "microbiome", "ggpubr", "viridis", "decontam", "gridExtra", "ggpubr", "lme4", "lmerTest", "writexl", "harrietr", "Maaslin2", "ggtext", "ggpmisc", "gridExtra", "gamm4", "reshape2", "kableExtra", "knitr")
YAML_header <-
'---
title: "Host-DNA depletion 1: data wrangling"
author: "Minsik Kim"
date: "2032.03.10"
output:
rmdformats::downcute:
downcute_theme: "chaos"
code_folding: hide
fig_width: 6
fig_height: 6
---'
seed <- "20230314"
#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
#Loads libraries, file paths, and other document options.
#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
FUN.LineZero.Boot <- function() {
.libPaths(path_library)
require(pacman)
pacman::p_load(c("knitr", "rmarkdown", "rmdformats", "yaml"))
knitr::opts_knit$set(root.dir = path_working)
str_libraries |> unique() |> sort() -> str_libraries
pacman::p_load(char = str_libraries)
set.seed(seed)
}
<<<<<<< HEAD
FUN.LineZero.Boot()
##
## The downloaded binary packages are in
## /var/folders/z2/kwm76s1n3jj2sjznjyj9thqc0000gn/T//Rtmp7e3TsQ/downloaded_packages
#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
=======
FUN.LineZero.Boot()
#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
>>>>>>> a74da024aec8928028c246dbcfe35b96e9691229
#Outputs R environment report.
#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
FUN.LineZero.Report <- function() {
cat("Line Zero Environment:\n\n")
paste("R:", pacman::p_version(), "\n") |> cat()
cat("Libraries:\n")
for (str_libraries in str_libraries) {
paste(
" ", str_libraries, ": ", pacman::p_version(package = str_libraries),
"\n", sep = ""
) |> cat()
}
paste("\nOperating System:", pacman::p_detectOS(), "\n") |> cat()
paste(" Library Path:", path_library, "\n") |> cat()
paste(" Working Path:", path_working, "\n") |> cat()
paste("Seed:", seed, "\n\n") |> cat()
cat("YAML Header:\n")
cat(YAML_header)
}
FUN.LineZero.Report()
## Line Zero Environment:
##
## R: 4.2.2
## Libraries:
## readxl: 1.4.2
<<<<<<< HEAD
## phyloseq: 1.42.0
=======
## phyloseq: 1.40.0
>>>>>>> a74da024aec8928028c246dbcfe35b96e9691229
## tidyverse: 2.0.0
## pacman: 0.5.1
## yaml: 2.3.7
## ggplot2: 3.4.1
## vegan: 2.6.4
<<<<<<< HEAD
## microbiome: 1.20.0
## ggpubr: 0.6.0
## viridis: 0.6.2
## decontam: 1.18.0
=======
## microbiome: 1.18.0
## ggpubr: 0.6.0
## viridis: 0.6.2
## decontam: 1.16.0
>>>>>>> a74da024aec8928028c246dbcfe35b96e9691229
## gridExtra: 2.3
## ggpubr: 0.6.0
## lme4: 1.1.31
## lmerTest: 3.1.3
## writexl: 1.4.2
## harrietr: 0.2.3
<<<<<<< HEAD
## Maaslin2: 1.12.0
=======
## Maaslin2: 1.10.0
>>>>>>> a74da024aec8928028c246dbcfe35b96e9691229
## ggtext: 0.1.2
## ggpmisc: 0.5.2
## gridExtra: 2.3
## gamm4: 0.2.6
## reshape2: 1.4.4
## kableExtra: 1.3.4
## knitr: 1.42
##
## Operating System: Darwin
## Library Path: /Library/Frameworks/R.framework/Resources/library
## Working Path: /Users/minsikkim/Dropbox (Partners HealthCare)/Project_SICAS2_microbiome/5_Scripts/MGK/Host_depletion_git
## Seed: 20230314
##
## YAML Header:
## ---
## title: "Host-DNA depletion 1: data wrangling"
## author: "Minsik Kim"
## date: "2032.03.10"
## output:
## rmdformats::downcute:
## downcute_theme: "chaos"
## code_folding: hide
## fig_width: 6
## fig_height: 6
## ---
#Script description#
*1. Loading data 1.1. phyloseq obejct 1.2. qPCR data (controls)
QC Q1. How many samples failed sequencing Q2. How were changes in read stats and host DNA proportion? Q3. How were the extraction controls Q4. Prevalence / abundance filtering - red flag
Exploratory analysis *
data input required:
Meta data
qPCR - bacteria
qPCR - human
qPCR host %
Raw reads
final reads
sequencing host %
library prep failure status
Raw reads
subject_id
treatment
sample_type
subject_id
Sequencing result
samples
controls
Analysis
calculation of alpha-diversity indices
qPCR, by smaple type and treatment A. Total DNA B. Host DNA C. Bacterial DNA D. Host %
Sequencing results (treatment, raw reads, host reads, % host, final reads) –> both stratified and nonstratified library prep % vs sample_type + treatment + sample_type * treatment + subject_id (binary) log10(Final reads) vs sample_type + treatment + sample_type * treatment + subject_id (linear or binary, see regular and cumulative distribution) Host DNA % vs sample_type + treatment + sample_type * treatment + subject_id (linear or binary, see regular and cumulative distribution) –> both stratified and nonstratified
LM of taxa alpha diversity (richness) –> both stratified and nonstratified
Permanova (Taxa dist ~ log10(final reads) + sample_type + treatment + sample_type * treatment + subject_id) –> both stratified and nonstratified
DA analysis for taxa, by sample type and treatment –> both stratified and nonstratified –> look at level groups as well - phylum, etc.
Decontam - stratified by treatment –> both stratified and nonstratified
LM of taxa alpha diversity (Shannon) –> both stratified and nonstratified
LM of taxa alpha diversity (InvSimpson) –> both stratified and nonstratified
LM of taxa alpha diversity (BPI) –> both stratified and nonstratified
LM of function alpha diversity (richness) –> both stratified and nonstratified
LM of function alpha diversity (Shannon) –> both stratified and nonstratified
LM of function alpha diversity (InvSimpson) –> both stratified and nonstratified
LM of function alpha diversity (BPI) –> both stratified and nonstratified
permanova of function alpha diversity –> both stratified and nonstratified
#Loading data#
# Loading files -----------------------------------------------------------
#loading tidy phyloseq object
phyloseq <- read_rds("/Users/minsikkim/Dropbox (Partners HealthCare)/Project_SICAS2_microbiome/4_Data/2_Tidy/Phyloseq/PHY_20221129_MGK_host_tidy_tax.rds")
####Need to be removed after adding sequenicing data
data_qPCR <- read_csv("/Users/minsikkim/Dropbox (Partners HealthCare)/Project_SICAS2_microbiome/7_Manuscripts/2022_MGK_Host_Depletion/Tables/DAT_20230122_MGK_host_control_qPCR.csv")
#qPCR data of all the samples sent out sequencing
data_qPCR <- subset(data_qPCR, data_qPCR$baylor_other_id %in% c(sample_names(phyloseq$phyloseq_count)) | data_qPCR$baylor_other_id %in% c(read_excel("/Users/minsikkim/Dropbox (Partners HealthCare)/Project_Baylor/3_Documentation/Communications/2023-01-24_baylor_shipping_sicas2_nasal_host_depleted/CMMR_MetadataCapture_20230124_LaiP-PQ00430_SICAS2_NS.xlsx", skip = 27) %>% data.frame %>% .$`Optional..............secondary.ID`))
data_qPCR <- subset(data_qPCR, data_qPCR$sample_type %in% c("BAL", "nasal_swab", "Sputum", "neg_depletion", "pos_depletion"))
data_qPCR$sample_type <- factor(data_qPCR$sample_type, levels = c("BAL", "nasal_swab", "Sputum", "pos_depletion", "neg_depletion"),
labels = c("BAL", "Nasal swab", "Sputum", "P depletion", "N depletion"))
data_qPCR$treatment <- factor(data_qPCR$treatment, levels = c("control", "lyPMA", "benzonase", "host_zero", "molysis", "qiaamp"),
labels = c("Control", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp"))
#sample data loading
sample_data <- sample_data(phyloseq$phyloseq_count)
- QC Q1. How many samples failed sequencing Q2. How were changes in read stats and host DNA proportion? Q3. How were the extraction controls Q4. Decontam - were there any contaminants?
Q1. How many samples failed in sequencing?
# Initail QC --------------------------------------------------------------
#Quesetions - QC
#Q0. How many samples failed in sequencing
## figures -----raw data
sample_data %>%
subset(., !is.na(.$subject_id)) %>%
data.frame() %>%
mutate(host_seq_percent = sequencing_host_prop * 100,
.after = sequencing_host_prop,) %>%
gather(feature, value, Raw_reads:host_seq_percent) %>%
group_by(feature, sample_type) %>%
subset(., .$feature %in% c("Raw_reads", "Host_mapped", "Final_reads", "host_seq_percent")) %>%
mutate(feature = factor(feature, levels = c("Raw_reads", "Host_mapped", "Final_reads", "host_seq_percent"), labels = c("Raw reads", "Host mapped", "Final reads", "Host %"))) %>%
ggplot(aes(x = value, fill = treatment)) +
geom_histogram(bins = 97) +
theme(legend.position = "top") +
guides(fill=guide_legend(title="Treatment", nrow = 1)) +
facet_grid(sample_type~feature, scales = "free") +
ggtitle("log10 transfromed histrogram")
## figures -----log10
sample_data %>%
subset(., !is.na(.$subject_id)) %>%
data.frame() %>%
mutate(host_seq_percent = sqrt(sequencing_host_prop),
.after = sequencing_host_prop,) %>%
gather(feature, value, Raw_reads:host_seq_percent) %>%
group_by(feature, sample_type) %>%
subset(., .$feature %in% c("Raw_reads", "Host_mapped", "Final_reads", "host_seq_percent")) %>%
mutate(feature = factor(feature, levels = c("Raw_reads", "Host_mapped", "Final_reads", "host_seq_percent"), labels = c("Raw reads", "Host mapped", "Final reads", "Host %"))) %>%
ggplot(aes(x = log10(value), fill = treatment)) +
geom_histogram(bins = 97) +
theme(legend.position = "top") +
guides(fill=guide_legend(title="Treatment", nrow = 1)) +
facet_grid(sample_type~feature, scales = "free") +
ggtitle("log10 for reads, sqrt for host % transfromed histrogram")
# no somaple showed 0 reads
ggarrange(ggplot(sample_data %>% subset(., !is.na(.$subject_id)) %>% data.frame(), aes(x = Final_reads, fill = treatment)) +
geom_histogram(bins = 97) +
facet_wrap(~sample_type) +
theme_classic(base_family = "serif") +
ggtitle("Histogram of final reads by sample type and treatment") +
scale_fill_manual(values = c("#e41a1c", "#377eb8", "#4daf4a", "#984ea3", "#ff7f00", "#ffff33"), name = "Treatment"),
ggplot(sample_data %>% subset(., !is.na(.$subject_id)) %>% data.frame(), aes(x = log10(Final_reads), fill = treatment)) +
geom_histogram(bins = 97) +
facet_wrap(~sample_type) +
theme_classic(base_family = "serif") +
ggtitle("Histogram of log10(final reads) by sample type and treatment") +
scale_fill_manual(values = c("#e41a1c", "#377eb8", "#4daf4a", "#984ea3", "#ff7f00", "#ffff33"), name = "Treatment"),
common.legend = T, ncol = 1)
hist((log10((phyloseq$phyloseq_count %>% otu_table %>% colSums()) + 1)),100, main = "Histogram of total reads (sum of OTU table)") # 2 samples showed 0 total reads (sum of otu_table)
#how were the samples failed in library prep?
sample_data %>% data.frame %>% mutate(total_read = phyloseq$phyloseq_count %>% otu_table %>% colSums()) %>% rownames_to_column(var = "baylor_other_id") %>%
ggplot(aes(x = reorder(baylor_other_id, -total_read),
y = log10(total_read + 1),
col = lib_failed)) +
geom_point() +
theme_classic(base_family = "serif") +
theme(axis.title.y = element_markdown(), axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size = 4)) +
ylab("log<sub>10</sub>(Sum of OTU table reads)") +
xlab("Sample ID") +
guides(col=guide_legend(title="Library failed")) +
ggtitle("Sum of OTU reads by library failure status")
sample_data %>% data.frame %>% filter(phyloseq$phyloseq_count %>% otu_table %>% colSums() == 0) # two BAL sampels showed 0 total reads
sample_data(phyloseq$phyloseq_count) %>% data.frame() %>% subset(., .$lib_failed)
Results:
2.1.1 Modeling final read should be conducted with both linear and logistic regression. Host % need no transformation, and should be tested for both linear and binomial methods
2.1.2 13 samples failed in library prep
2.1.3. Two BAL sampels showed 0 total reads
2.1.4. Sequencing fail ≠ library prep fail
Comments from Baylor:
Q: What was the lab’s criteria for deciding which samples failed library prep.? There were 13 samples that you pointed as failed but their sequencing result actually looks pretty good (ie similar to samples that didn’t fail library prep)
A: To determine whether a library attempt “passed or failed” the lab looks at the picogreen concentrations and a library fragment size distribution trace. The trace files are an output from either the Fragment Analyzer or TapeStation (a copy of the trace files for PQ00331 is attached). If a sample has a background level pico concentration and no discernable fragment concentration on the trace (i.e. a flat trace line) it is considered failed library. If the sample is below the level of detection of our standard library QC methods, it is considered failure. It’s still possible that there is some small amounts of library in those samples that were successfully sequenced, but often those samples do not generate a meaningful amount of sequence data.
QC2 table
#sequencing result by sample type and control (1/0)
options(dplyr.summarise.inform = FALSE)
sample_data %>% data.frame() %>%
group_by(sample_type, treated) %>%
summarise(N = n(),
`Raw reads<br>(median [IQR])<br>[reads x 10<sup>7</sup>]` = paste(format(round(median(Raw_reads/10000000),2),nsmall = 2, big.mark = ","), " [", format(round(quantile(Raw_reads/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Raw_reads/10000000, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
`Host reads<br>(median [IQR])<br>[reads x 10<sup>7</sup>]` = paste(format(round(median(Host_mapped/10000000),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(Host_mapped/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Host_mapped/10000000, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
`Host reads proportion<br>(median [IQR])<br>[%]` = paste(format(round(median(sequencing_host_prop * 100),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(sequencing_host_prop * 100, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(sequencing_host_prop * 100, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
`Final reads<br>(median [IQR])<br>[reads x 10<sup>7</sup>]` = paste(format(round(median(Final_reads/10000000),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(Final_reads/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Final_reads/10000000, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
) %>%
rename(`Sample type` = sample_type, Treated = treated) %>%
data.frame(check.names = F) %>% mutate_all(linebreak) %>% kbl(format = "html", escape = F) %>% kable_styling(full_width = 0, html_font = "serif")
| Sample type | Treated | N |
Raw reads (median [IQR]) [reads x 107] |
Host reads (median [IQR]) [reads x 107] |
Host reads proportion (median [IQR]) [%] |
Final reads (median [IQR]) [reads x 107] |
|---|---|---|---|---|---|---|
| BAL | 0 | 5 | 15.73 [6.35, 15.92] | 12.92 [5.21, 12.94] | 82.14 [82.12, 82.18] | 0.03 [0.03, 0.04] |
| BAL | 1 | 25 | 6.17 [4.57, 17.43] | 4.65 [2.78, 12.80] | 75.59 [69.38, 78.44] | 0.17 [0.10, 0.37] |
| Nasal swab | 0 | 10 | 13.09 [7.73, 16.93] | 10.05 [6.11, 13.04] | 77.34 [76.89, 79.45] | 0.48 [0.10, 0.87] |
| Nasal swab | 1 | 25 | 4.08 [0.99, 6.40] | 0.81 [0.26, 1.36] | 27.46 [13.51, 64.58] | 0.97 [0.17, 3.42] |
| Sputum | 0 | 5 | 8.59 [8.25, 9.27] | 6.87 [6.69, 7.50] | 80.61 [79.92, 80.89] | 0.06 [0.06, 0.09] |
| Sputum | 1 | 25 | 12.23 [10.34, 13.73] | 7.71 [3.76, 8.82] | 58.64 [39.71, 74.93] | 1.16 [0.47, 4.19] |
| Neg. ext. | 0 | 1 | 0.02 [0.02, 0.02] | 0.00 [0.00, 0.00] | 3.13 [3.13, 3.13] | 0.02 [0.02, 0.02] |
| Pos. ext. | 0 | 1 | 11.87 [11.87, 11.87] | 0.00 [0.00, 0.00] | 0.02 [0.02, 0.02] | 9.96 [9.96, 9.96] |
sample_data %>% data.frame() %>%
#dplyr::filter(sample_type %in% c("Sputum", "nasal_swab", "BAL")) %>%
group_by (sample_type, treatment) %>%
summarise(N = n(),
`Raw reads<br>(median [IQR])<br>[reads x 10<sup>7</sup>]` = paste(format(round(median(Raw_reads/10000000),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(Raw_reads/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Raw_reads/10000000, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
`Host reads<br>(median [IQR])<br>[reads x 10<sup>7</sup>]` = paste(format(round(median(Host_mapped/10000000),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(Host_mapped/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Host_mapped/10000000, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
`Host reads proportion<br>(median [IQR])<br>[%]` = paste(format(round(median(sequencing_host_prop * 100),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(sequencing_host_prop * 100, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(sequencing_host_prop * 100, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
`Final reads<br>(median [IQR])<br>[reads x 10<sup>7</sup>]` = paste(format(round(median(Final_reads/10000000),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(Final_reads/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Final_reads/10000000, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
) %>% data.frame(check.names = F) %>%
arrange(sample_type, treatment) %>%
rename(`Sample type` = sample_type, Treatment = treatment) %>%
mutate_all(linebreak) %>% kbl(format = "html", escape = F) %>% kable_styling(full_width = 0, html_font = "serif")
| Sample type | Treatment | N |
Raw reads (median [IQR]) [reads x 107] |
Host reads (median [IQR]) [reads x 107] |
Host reads proportion (median [IQR]) [%] |
Final reads (median [IQR]) [reads x 107] |
|---|---|---|---|---|---|---|
| BAL | Control | 5 | 15.73 [6.35, 15.92] | 12.92 [5.21, 12.94] | 82.14 [82.12, 82.18] | 0.03 [0.03, 0.04] |
| BAL | lyPMA | 5 | 5.72 [3.59, 13.41] | 4.65 [2.79, 10.90] | 80.73 [77.85, 81.18] | 0.06 [0.04, 0.10] |
| BAL | Benzonase | 5 | 18.59 [16.20, 23.63] | 14.77 [12.80, 18.16] | 79.02 [77.89, 79.45] | 0.17 [0.16, 0.22] |
| BAL | Host zero | 5 | 4.57 [2.32, 4.71] | 2.69 [1.61, 2.93] | 66.98 [57.17, 68.95] | 0.24 [0.13, 0.82] |
| BAL | Molysis | 5 | 4.76 [3.57, 4.86] | 2.78 [1.39, 3.61] | 74.39 [74.29, 75.59] | 0.29 [0.13, 1.56] |
| BAL | QIAamp | 5 | 17.19 [15.35, 17.43] | 11.87 [10.79, 12.22] | 74.88 [71.08, 77.31] | 0.26 [0.10, 1.02] |
| Nasal swab | Control | 10 | 13.09 [7.73, 16.93] | 10.05 [6.11, 13.04] | 77.34 [76.89, 79.45] | 0.48 [0.10, 0.87] |
| Nasal swab | lyPMA | 5 | 0.98 [0.85, 1.24] | 0.63 [0.28, 0.88] | 71.37 [28.73, 74.68] | 0.07 [0.06, 0.08] |
| Nasal swab | Benzonase | 5 | 5.75 [4.95, 6.57] | 3.66 [1.29, 5.05] | 64.58 [63.71, 76.83] | 0.28 [0.26, 1.04] |
| Nasal swab | Host zero | 5 | 2.83 [1.42, 6.42] | 0.49 [0.03, 0.81] | 7.67 [2.33, 25.73] | 2.43 [0.97, 5.03] |
| Nasal swab | Molysis | 5 | 0.99 [0.63, 4.08] | 0.42 [0.06, 0.64] | 41.04 [4.31, 64.24] | 0.32 [0.17, 2.53] |
| Nasal swab | QIAamp | 5 | 6.40 [6.40, 6.80] | 0.86 [0.86, 1.17] | 17.24 [13.51, 19.57] | 4.63 [4.50, 4.67] |
| Sputum | Control | 5 | 8.59 [8.25, 9.27] | 6.87 [6.69, 7.50] | 80.61 [79.92, 80.89] | 0.06 [0.06, 0.09] |
| Sputum | lyPMA | 5 | 10.98 [5.22, 12.78] | 8.82 [3.76, 10.44] | 76.33 [72.02, 80.29] | 0.25 [0.15, 0.44] |
| Sputum | Benzonase | 5 | 10.76 [10.34, 10.82] | 7.81 [7.75, 8.24] | 74.93 [72.13, 75.33] | 0.47 [0.45, 0.59] |
| Sputum | Host zero | 5 | 13.14 [7.64, 13.95] | 4.39 [3.80, 7.71] | 49.71 [31.45, 55.49] | 2.91 [2.36, 3.67] |
| Sputum | Molysis | 5 | 12.59 [10.84, 13.73] | 2.98 [1.83, 4.28] | 27.49 [14.62, 28.19] | 6.11 [5.56, 8.37] |
| Sputum | QIAamp | 5 | 12.35 [12.23, 12.85] | 9.08 [8.41, 9.27] | 71.76 [56.69, 73.50] | 1.16 [1.13, 3.89] |
| Neg. ext. | Control | 1 | 0.02 [0.02, 0.02] | 0.00 [0.00, 0.00] | 3.13 [3.13, 3.13] | 0.02 [0.02, 0.02] |
| Pos. ext. | Control | 1 | 11.87 [11.87, 11.87] | 0.00 [0.00, 0.00] | 0.02 [0.02, 0.02] | 9.96 [9.96, 9.96] |
QC2 figure
# Summary figures - facet and z-score -------------------------------------
sample_data %>%
subset(., !is.na(.$subject_id)) %>%
data.frame() %>%
gather(feature, value, Raw_reads:sequencing_host_prop) %>%
group_by(feature, sample_type) %>%
subset(., .$feature %in% c("Raw_reads", "Host_mapped", "Final_reads", "sequencing_host_prop")) %>%
mutate(z_score = scale(value),
feature = factor(feature, levels = c("Raw_reads", "Host_mapped", "Final_reads", "sequencing_host_prop"), labels = c("Raw reads", "Host mapped", "Final reads", "Host %"))) %>%
ggplot(aes(x = treatment, y = z_score, fill = treatment)) +
geom_boxplot(lwd = 0.2) +
guides(fill=guide_legend(title="Treatment", nrow = 1)) +
facet_grid(sample_type~feature) +
xlab("Treatment") +
ylab("Z score") +
theme_classic(base_family = "serif", base_size = 14) +
guides( x = guide_axis(angle = 90)) +
theme(legend.position = "top") +
scale_fill_manual(values = c("#e41a1c", "#377eb8", "#4daf4a", "#984ea3", "#ff7f00", "#ffff33"), name = "Treatment") #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
Results:
2.2.1. There were no differences between raw reads, however, final reads increased after some treatment, and host DNA proportion decreased
2.2.2. The changes in final reads were varied by treatment and sample type
comment: Should I normalize z-score to mean value of control?
Q3. How was the extraction controls, both pos and neg?
#Loading theoretical mock community
zymo_mock <- read_excel("/Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_sicas2/data_raw/DAR_20210929_zymo_mock_data.xlsx") %>%
data.frame(row.names = T) %>% rename(mock_theoretical = Mock) %>% mutate(mock_theoretical = mock_theoretical/100) %>%
merge_phyloseq(otu_table(., taxa_are_rows = T), tax_table(phyloseq$phyloseq_count))
phyloseq_control_rel <- rbind(c("mock_theoretical", "mock")) %>% data.frame() %>%
column_to_rownames(var = "X1") %>% rename(sample_type = X2) %>% #making sample_data of mock community
merge_phyloseq(sample_data(.), zymo_mock) %>%
merge_phyloseq(., subset_samples(phyloseq$phyloseq_rel, sample_type == "Pos. ext." | sample_type == "Neg. ext." )) #adding data of controls
#Species richness of each control groups
cat("Species richness\n")
## Species richness
rowSums(t(otu_table(phyloseq_control_rel)) != 0)
## mock_theoretical 20220606_Pos 20220606_Neg
## 10 41 3
sample_data(phyloseq_control_rel)$sample_type <- factor(sample_data(phyloseq_control_rel)$sample_type,
levels = c("mock", "2", "1"),
labels = c("Mock, 10 spp.", "Pos., 41 spp.", "Neg., 3 spp."))
#Manipulating phyloseq - only top 10
tax_table(phyloseq_control_rel) %>%
cbind(species10 = "[Others]") %>%
{top10species <- head(taxa_sums(phyloseq_control_rel), 10) %>% names()
.[top10species, "species10"] <- as.character(.[top10species, "Species"])
.[, 8] <- .[, 8] %>% gsub("s__", "", .) %>% gsub("_", " ", .) %>% paste("<i>", ., "</i>", sep = "")
phyloseq_temp <- phyloseq_control_rel
tax_table(phyloseq_temp) <- tax_table(.)
phyloseq_temp
} %>%
plot_bar(., fill="species10") +
ylab("Relative abundancne") +
theme_classic(base_size = 11, base_family = "serif") +
ggtitle("Species richness of control data") +
theme(legend.text = element_markdown()) +
guides(fill=guide_legend(title="Top 10 species")) +
facet_wrap (~ sample_type, scales= "free_x", nrow=1)
#there could be opportunistic pathogens...
Results
2.3.1. Negative control showed minimal number of possible contaminants
2.3.2. Positive control contained various contaminants
Results:
2.2.1. There were no differences between raw reads, however, final reads increased after some treatment, and host DNA proportion decreased
2.2.2. The changes in final reads were varied by treatment and sample type
Q4. Prevalence filtering - red flag
#Calculation of sample prevalence, standard deviation, median abundance across all samples for all bugs and making into a table.
#
#• In initial analysis we will not perform prevalence or abundance filtering (though we may consider this for secondary differential abundance analyses to manage p (features) > n (sample size) problem and issues with multiple hypothesis correction)
taxa_qc <- data.frame("species" = otu_table(phyloseq$phyloseq_rel) %>% t() %>% colnames(),
"prevalence" = ifelse(phyloseq$phyloseq_count %>% otu_table() > 0, 1, 0) %>% t() %>% colSums(), #Prevalence of taxa
"mean_rel_abd" = phyloseq$phyloseq_rel %>% otu_table() %>% t() %>% colMeans(na.rm = T) #mean relativ abundacne
)
hist(log10(taxa_qc$prevalence), xlab = "log10(Taxa prevalence)", main = "Histogram of prevalence of taxa")
hist(log10(taxa_qc$mean_rel_abd), xlab = "log10(Mean relative abundance)", main = "Histogram of mean relative abundance")
• In initial analysis we will not perform prevalence or abundance filtering (though we may consider this for secondary differentialabundance analyses to manage p (features) > n (sample size) problem and issues with multiple hypothesis correction) • Calculation of sample prevalence, standard deviation, median abundance across all samples for all bugs and making into a table.
Although we don’t consider the prevalence of abundance at this time, we can consider their red-flags after running the DA analysis
- red flag
red_flag_taxa <- data.frame(species = taxa_qc$species,
red_flag_prev_abd = ifelse(taxa_qc$prevalence < otu_table(phyloseq$phyloseq_rel) %>% t %>% rownames() %>% length * 0.05 & taxa_qc$mean_rel_abd < quantile(taxa_qc$mean_rel_abd, 0.75), 1, 0)
)
#Analysis#
A0. calculation of alpha-diversity indices
alpha_diversity <- function(data) {
otu_table <- otu_table(data) %>% .[colSums(.) !=0]
S.obs <- rowSums(t(otu_table) != 0)
sample_data <- sample_data(data)
data_evenness <- vegan::diversity(t(otu_table)) / log(vegan::specnumber(t(otu_table))) # calculate evenness index using vegan package
data_shannon <- vegan::diversity(t(otu_table), index = "shannon") # calculate Shannon index using vegan package
data_hill <- exp(data_shannon) # calculate Hills index
data_dominance <- microbiome::dominance(otu_table, index = "all", rank = 1, aggregate = TRUE) # dominance (Berger-Parker index), etc.
data_invsimpson <- vegan::diversity(t(otu_table), index = "invsimpson") # calculate Shannon index using vegan package
alpha_diversity <- cbind(S.obs, data_shannon, data_hill, data_invsimpson, data_evenness,data_dominance) # combine all indices in one data table
sample_data <- merge(data.frame(sample_data), alpha_diversity, by = 0, all = T) %>% column_to_rownames(var = "Row.names")
}
sample_data(phyloseq$phyloseq_rel) <- sample_data(alpha_diversity(phyloseq$phyloseq_count))
sample_data(phyloseq$phyloseq_count) <- sample_data(alpha_diversity(phyloseq$phyloseq_count))
sample_data(phyloseq$phyloseq_path_rpkm) <- sample_data(alpha_diversity(phyloseq$phyloseq_path_rpkm))
A1. qPCR, by smaple type and treatment
A. Total DNA
B. Host DNA
C. Bacterial DNA
D. Host %
#2A: Change in total DNA (qPCR)
f2a <- ggplot(data_qPCR, aes(x = sample_type, y = log10(DNA_host_nondil + DNA_bac_nondil))) +
geom_boxplot(aes(fill = treatment), lwd = 0.2) +
scale_fill_manual(values = c("#e41a1c", "#377eb8", "#4daf4a", "#984ea3", "#ff7f00", "#ffff33")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
ylab("log<sub>10</sub>(qPCR total DNA)<br>(ng/μL)") +
xlab("Sample type") +
theme_classic (base_size = 12, base_family = "serif") +
labs(tag = "A") +
scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
theme(plot.tag = element_text(size = 15), axis.title.y = element_markdown()) + # Plot title size
guides(fill = guide_legend(nrow = 1, title = "Treatment"))
#2B: Change in human DNA (qPCR)
f2b <- ggplot(data_qPCR, aes(x = sample_type, y = log10(DNA_host_nondil))) +
geom_boxplot(aes(fill = treatment), lwd = 0.2) +
scale_fill_manual(values = c("#e41a1c", "#377eb8", "#4daf4a", "#984ea3", "#ff7f00", "#ffff33")) +
ylab("log<sub>10</sub>(qPCR host DNA)<br>(ng/μL)") +
xlab("Sample type") +
theme_classic (base_size = 12, base_family = "serif")+
labs(tag = "B") +
scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
theme(plot.tag = element_text(size = 15), axis.title.y = element_markdown()) + # Plot title size
guides(fill = guide_legend(nrow = 1, title = "Treatment"))
#2C: Change in 16S DNA (qPCR)
f2c <- ggplot(data_qPCR, aes(x = sample_type, y = log10(DNA_bac_nondil))) +
geom_boxplot(aes(fill = treatment), lwd = 0.2) +
scale_fill_manual(values = c("#e41a1c", "#377eb8", "#4daf4a", "#984ea3", "#ff7f00", "#ffff33")) +
ylab("log<sub>10</sub>(qPCR bacterial DNA)<br>(ng/μL)") +
xlab("Sample type") +
theme_classic (base_size = 12, base_family = "serif")+
labs(tag = "C") +
scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
theme(plot.tag = element_text(size = 15), axis.title.y = element_markdown()) + # Plot title size
guides(fill = guide_legend(nrow = 1, title = "Treatment"))
#2D. Change in % host (qPCR)
f2d <- ggplot(data_qPCR, aes(x = sample_type, y = log10(host_proportion * 100))) +
geom_boxplot(aes(fill = treatment), lwd = 0.2) +
scale_fill_manual(values = c("#e41a1c", "#377eb8", "#4daf4a", "#984ea3", "#ff7f00", "#ffff33")) +
ylab("log<sub>10</sub>(host DNA) (%)") +
xlab("Sample type") +
theme_classic (base_size = 12, base_family = "serif") +
labs(tag = "D") +
scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
theme(plot.tag = element_text(size = 15), axis.title.y = element_markdown()) + # Plot title size
guides(fill = guide_legend(nrow = 1, title = "Treatment"))
#output for markdown
ggarrange(f2a, f2b, f2c, f2d, common.legend = T , align = "hv")
A2. Sequencing results (treatment, raw reads, host reads, % host, final reads)
--> both stratified and nonstratified
library prep % vs sample_type + treatment + sample_type * treatment + subject_id (binary)
log10(Final reads) vs sample_type + treatment + sample_type * treatment + subject_id (linear or binary, see regular and cumulative distribution)
Host DNA % vs sample_type + treatment + sample_type * treatment + subject_id (linear or binary, see regular and cumulative distribution)
--> both stratified and nonstratified
cat("Binomial - library failed? \n\n")
## Binomial - library failed?
glmer(lib_failed ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
summary %>% .$coefficients %>% data.frame(check.names = F) %>%
mutate(` ` = case_when(abs(`t value`) > 2 ~ "*", .default = " ")) %>%
kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")
| Estimate | Std. Error | t value | ||
|---|---|---|---|---|
| (Intercept) | 0.0000000 | 0.1167011 | 0.0000000 | |
| sample_typeNasal swab | -0.0204626 | 0.1470919 | -0.1391145 | |
| sample_typeSputum | 0.0000000 | 0.1650403 | 0.0000000 | |
| treatmentlyPMA | 0.2000000 | 0.1619675 | 1.2348160 | |
| treatmentBenzonase | 0.0000000 | 0.1619675 | 0.0000000 | |
| treatmentHost zero | 0.2000000 | 0.1619675 | 1.2348160 | |
| treatmentMolysis | 0.2000000 | 0.1619675 | 1.2348160 | |
| treatmentQIAamp | 0.0000000 | 0.1619675 | 0.0000000 | |
| sample_typeNasal swab:treatmentlyPMA | 0.5978181 | 0.2143680 | 2.7887474 |
|
| sample_typeSputum:treatmentlyPMA | -0.2000000 | 0.2290566 | -0.8731468 | |
| sample_typeNasal swab:treatmentBenzonase | -0.0021819 | 0.2143680 | -0.0101781 | |
| sample_typeSputum:treatmentBenzonase | 0.0000000 | 0.2290566 | 0.0000000 | |
| sample_typeNasal swab:treatmentHost zero | 0.1978181 | 0.2143680 | 0.9227971 | |
| sample_typeSputum:treatmentHost zero | -0.2000000 | 0.2290566 | -0.8731468 | |
| sample_typeNasal swab:treatmentMolysis | 0.6021819 | 0.2143680 | 2.8091035 |
|
| sample_typeSputum:treatmentMolysis | -0.2000000 | 0.2290566 | -0.8731468 | |
| sample_typeNasal swab:treatmentQIAamp | 0.0021819 | 0.2143680 | 0.0101781 | |
| sample_typeSputum:treatmentQIAamp | 0.0000000 | 0.2290566 | 0.0000000 |
cat("Binomial - Final reads? \n\n")
## Binomial - Final reads?
glmer(log10(Final_reads) ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
summary %>% .$coefficients %>% data.frame(check.names = F) %>%
mutate(` ` = case_when(abs(`t value`) > 2 ~ "*", .default = " ")) %>%
kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")
| Estimate | Std. Error | t value | ||
|---|---|---|---|---|
| (Intercept) | 5.4933063 | 0.2095210 | 26.2184086 |
|
| sample_typeNasal swab | 1.1246781 | 0.2789651 | 4.0316080 |
|
| sample_typeSputum | 0.3332900 | 0.2963074 | 1.1248117 | |
| treatmentlyPMA | 0.3521385 | 0.2733402 | 1.2882794 | |
| treatmentBenzonase | 0.8109049 | 0.2733402 | 2.9666509 |
|
| treatmentHost zero | 0.9501904 | 0.2733402 | 3.4762193 |
|
| treatmentMolysis | 1.0358039 | 0.2733402 | 3.7894313 |
|
| treatmentQIAamp | 1.0480985 | 0.2733402 | 3.8344107 |
|
| sample_typeNasal swab:treatmentlyPMA | -0.8774152 | 0.3621899 | -2.4225278 |
|
| sample_typeSputum:treatmentlyPMA | 0.1879364 | 0.3865614 | 0.4861747 | |
| sample_typeNasal swab:treatmentBenzonase | -0.6887274 | 0.3621899 | -1.9015642 | |
| sample_typeSputum:treatmentBenzonase | 0.0349655 | 0.3865614 | 0.0904526 | |
| sample_typeNasal swab:treatmentHost zero | -0.1157121 | 0.3621899 | -0.3194791 | |
| sample_typeSputum:treatmentHost zero | 0.7231403 | 0.3865614 | 1.8706997 | |
| sample_typeNasal swab:treatmentMolysis | -0.8411412 | 0.3621899 | -2.3223760 |
|
| sample_typeSputum:treatmentMolysis | 0.9543008 | 0.3865614 | 2.4686913 |
|
| sample_typeNasal swab:treatmentQIAamp | 0.0236242 | 0.3621899 | 0.0652260 | |
| sample_typeSputum:treatmentQIAamp | 0.3676901 | 0.3865614 | 0.9511817 |
cat("Linear - Final reads? \n\n")
## Linear - Final reads?
lmer(log10(Final_reads) ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
summary %>% .$coefficients %>% data.frame(check.names = F) %>%
mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*", .default = " ")) %>%
kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")
| Estimate | Std. Error | df | t value | Pr(>|t|) | ||
|---|---|---|---|---|---|---|
| (Intercept) | 5.4933063 | 0.2095210 | 62.79820 | 26.2184086 | 0.0000000 |
|
| sample_typeNasal swab | 1.1246781 | 0.2789651 | 37.26869 | 4.0316080 | 0.0002633 |
|
| sample_typeSputum | 0.3332900 | 0.2963074 | 62.79820 | 1.1248117 | 0.2649498 | |
| treatmentlyPMA | 0.3521385 | 0.2733402 | 67.65016 | 1.2882794 | 0.2020376 | |
| treatmentBenzonase | 0.8109049 | 0.2733402 | 67.65016 | 2.9666509 | 0.0041580 |
|
| treatmentHost zero | 0.9501904 | 0.2733402 | 67.65016 | 3.4762193 | 0.0008929 |
|
| treatmentMolysis | 1.0358039 | 0.2733402 | 67.65016 | 3.7894313 | 0.0003238 |
|
| treatmentQIAamp | 1.0480985 | 0.2733402 | 67.65016 | 3.8344107 | 0.0002788 |
|
| sample_typeNasal swab:treatmentlyPMA | -0.8774152 | 0.3621899 | 68.18953 | -2.4225278 | 0.0180749 |
|
| sample_typeSputum:treatmentlyPMA | 0.1879364 | 0.3865614 | 67.65016 | 0.4861747 | 0.6284145 | |
| sample_typeNasal swab:treatmentBenzonase | -0.6887274 | 0.3621899 | 68.18953 | -1.9015642 | 0.0614539 | |
| sample_typeSputum:treatmentBenzonase | 0.0349655 | 0.3865614 | 67.65016 | 0.0904526 | 0.9281949 | |
| sample_typeNasal swab:treatmentHost zero | -0.1157121 | 0.3621899 | 68.18953 | -0.3194791 | 0.7503399 | |
| sample_typeSputum:treatmentHost zero | 0.7231403 | 0.3865614 | 67.65016 | 1.8706997 | 0.0657128 | |
| sample_typeNasal swab:treatmentMolysis | -0.8411412 | 0.3621899 | 68.18953 | -2.3223760 | 0.0232052 |
|
| sample_typeSputum:treatmentMolysis | 0.9543008 | 0.3865614 | 67.65016 | 2.4686913 | 0.0160928 |
|
| sample_typeNasal swab:treatmentQIAamp | 0.0236242 | 0.3621899 | 68.18953 | 0.0652260 | 0.9481849 | |
| sample_typeSputum:treatmentQIAamp | 0.3676901 | 0.3865614 | 67.65016 | 0.9511817 | 0.3448982 |
cat("Binomial - Host DNA proportion? \n\n")
## Binomial - Host DNA proportion?
glmer(sequencing_host_prop ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
summary %>% .$coefficients %>% data.frame(check.names = F) %>%
mutate(` ` = case_when(abs(`t value`) > 2 ~ "*", .default = " ")) %>%
kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")
| Estimate | Std. Error | t value | ||
|---|---|---|---|---|
| (Intercept) | 0.8200713 | 0.0646272 | 12.6892671 |
|
| sample_typeNasal swab | -0.0378253 | 0.0799149 | -0.4733196 | |
| sample_typeSputum | -0.0293956 | 0.0913966 | -0.3216265 | |
| treatmentlyPMA | -0.0377605 | 0.0908962 | -0.4154243 | |
| treatmentBenzonase | -0.0324363 | 0.0908962 | -0.3568498 | |
| treatmentHost zero | -0.2000025 | 0.0908962 | -2.2003393 |
|
| treatmentMolysis | -0.1577230 | 0.0908962 | -1.7351990 | |
| treatmentQIAamp | -0.0928027 | 0.0908962 | -1.0209737 | |
| sample_typeNasal swab:treatmentlyPMA | -0.1931415 | 0.1202628 | -1.6059959 | |
| sample_typeSputum:treatmentlyPMA | 0.0106590 | 0.1285467 | 0.0829196 | |
| sample_typeNasal swab:treatmentBenzonase | -0.1301969 | 0.1202628 | -1.0826036 | |
| sample_typeSputum:treatmentBenzonase | -0.0432386 | 0.1285467 | -0.3363648 | |
| sample_typeNasal swab:treatmentHost zero | -0.3924910 | 0.1202628 | -3.2636117 |
|
| sample_typeSputum:treatmentHost zero | -0.1532698 | 0.1285467 | -1.1923278 | |
| sample_typeNasal swab:treatmentMolysis | -0.2637315 | 0.1202628 | -2.1929608 |
|
| sample_typeSputum:treatmentMolysis | -0.3863225 | 0.1285467 | -3.0053090 |
|
| sample_typeNasal swab:treatmentQIAamp | -0.5240678 | 0.1202628 | -4.3576894 |
|
| sample_typeSputum:treatmentQIAamp | -0.0370057 | 0.1285467 | -0.2878775 |
cat("Linear - Host DNA proportion? \n\n")
## Linear - Host DNA proportion?
lmer(sequencing_host_prop ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
summary %>% .$coefficients %>% data.frame(check.names = F) %>%
mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*", .default = " ")) %>%
kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")
| Estimate | Std. Error | df | t value | Pr(>|t|) | ||
|---|---|---|---|---|---|---|
| (Intercept) | 0.8200713 | 0.0646272 | 75.67093 | 12.6892671 | 0.0000000 |
|
| sample_typeNasal swab | -0.0378253 | 0.0799149 | 41.00023 | -0.4733196 | 0.6384957 | |
| sample_typeSputum | -0.0293956 | 0.0913966 | 75.67093 | -0.3216265 | 0.7486224 | |
| treatmentlyPMA | -0.0377605 | 0.0908962 | 68.23467 | -0.4154243 | 0.6791348 | |
| treatmentBenzonase | -0.0324363 | 0.0908962 | 68.23467 | -0.3568498 | 0.7223055 | |
| treatmentHost zero | -0.2000025 | 0.0908962 | 68.23467 | -2.2003393 | 0.0311715 |
|
| treatmentMolysis | -0.1577230 | 0.0908962 | 68.23467 | -1.7351990 | 0.0872202 | |
| treatmentQIAamp | -0.0928027 | 0.0908962 | 68.23467 | -1.0209737 | 0.3108732 | |
| sample_typeNasal swab:treatmentlyPMA | -0.1931415 | 0.1202628 | 68.80454 | -1.6059959 | 0.1128550 | |
| sample_typeSputum:treatmentlyPMA | 0.0106590 | 0.1285467 | 68.23467 | 0.0829196 | 0.9341582 | |
| sample_typeNasal swab:treatmentBenzonase | -0.1301969 | 0.1202628 | 68.80454 | -1.0826036 | 0.2827639 | |
| sample_typeSputum:treatmentBenzonase | -0.0432386 | 0.1285467 | 68.23467 | -0.3363648 | 0.7376281 | |
| sample_typeNasal swab:treatmentHost zero | -0.3924910 | 0.1202628 | 68.80454 | -3.2636117 | 0.0017148 |
|
| sample_typeSputum:treatmentHost zero | -0.1532698 | 0.1285467 | 68.23467 | -1.1923278 | 0.2372628 | |
| sample_typeNasal swab:treatmentMolysis | -0.2637315 | 0.1202628 | 68.80454 | -2.1929608 | 0.0316937 |
|
| sample_typeSputum:treatmentMolysis | -0.3863225 | 0.1285467 | 68.23467 | -3.0053090 | 0.0037094 |
|
| sample_typeNasal swab:treatmentQIAamp | -0.5240678 | 0.1202628 | 68.80454 | -4.3576894 | 0.0000449 |
|
| sample_typeSputum:treatmentQIAamp | -0.0370057 | 0.1285467 | 68.23467 | -0.2878775 | 0.7743131 |
A3. LM of taxa alpha diversity
--> both stratified and nonstratified
sample_data <- sample_data(phyloseq$phyloseq_count) %>% data.frame(check.names = F) %>% subset(., !is.nan(.$simpson))
cat("Species richness of all samples\n\n")
## Species richness of all samples
cat("S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample)\n\n")
## S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample)
lmer(S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|subject_id), data = sample_data) %>%
summary() %>%
.$coefficients %>%
data.frame(check.names = F) %>%
mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Estimate | Std. Error | df | t value | Pr(>|t|) | ||
|---|---|---|---|---|---|---|
| (Intercept) | -62.286176 | 18.946348 | 73.91996 | -3.2875031 | 0.0015486 |
|
| sample_typeNasal swab | -9.894010 | 13.143390 | 19.76996 | -0.7527746 | 0.4604558 | |
| sample_typeSputum | 5.066489 | 11.833705 | 23.07241 | 0.4281405 | 0.6725190 | |
| treatmentlyPMA | -3.032607 | 7.992068 | 64.29757 | -0.3794521 | 0.7056025 | |
| treatmentBenzonase | -4.719045 | 8.051585 | 65.30664 | -0.5861014 | 0.5598277 | |
| treatmentHost zero | -2.064590 | 8.208256 | 65.44300 | -0.2515260 | 0.8021956 | |
| treatmentMolysis | 10.862491 | 8.314322 | 65.52729 | 1.3064795 | 0.1959566 | |
| treatmentQIAamp | -2.491588 | 8.330143 | 65.53938 | -0.2991050 | 0.7658063 | |
| log10(Final_reads) | 12.532134 | 3.113656 | 67.56866 | 4.0248934 | 0.0001465 |
|
| sample_typeNasal swab:treatmentlyPMA | 4.202591 | 10.419319 | 64.56848 | 0.4033461 | 0.6880263 | |
| sample_typeSputum:treatmentlyPMA | 34.664316 | 10.599221 | 64.21630 | 3.2704588 | 0.0017287 |
|
| sample_typeNasal swab:treatmentBenzonase | 1.775047 | 10.042806 | 64.92798 | 0.1767481 | 0.8602564 | |
| sample_typeSputum:treatmentBenzonase | 62.518484 | 10.357810 | 64.45145 | 6.0358784 | 0.0000001 |
|
| sample_typeNasal swab:treatmentHost zero | 4.793943 | 9.787379 | 64.68743 | 0.4898086 | 0.6259264 | |
| sample_typeSputum:treatmentHost zero | 92.094185 | 10.559614 | 64.43191 | 8.7213589 | 0.0000000 |
|
| sample_typeNasal swab:treatmentMolysis | -7.289176 | 10.178107 | 65.15977 | -0.7161623 | 0.4764500 | |
| sample_typeSputum:treatmentMolysis | 87.197251 | 10.723044 | 64.48528 | 8.1317631 | 0.0000000 |
|
| sample_typeNasal swab:treatmentQIAamp | -3.726532 | 9.774308 | 64.66718 | -0.3812579 | 0.7042616 | |
| sample_typeSputum:treatmentQIAamp | 70.548734 | 10.400892 | 64.40670 | 6.7829501 | 0.0000000 |
|
#Species richness - stratified
cat("Species richness at NS\n\n")
## Species richness at NS
cat("S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(S.obs ~ treatment + log10 (Final_reads) + (1|subject_id), data = subset(sample_data, sample_data$sample_type == "Nasal swab")) %>%
summary() %>%
.$coefficients %>%
data.frame(check.names = F) %>%
mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Estimate | Std. Error | df | t value | Pr(>|t|) | ||
|---|---|---|---|---|---|---|
| (Intercept) | -26.0233860 | 9.340261 | 28 | -2.7861520 | 0.0094654 |
|
| treatmentlyPMA | -2.1153025 | 2.324810 | 28 | -0.9098818 | 0.3706512 | |
| treatmentBenzonase | -1.7637014 | 2.206148 | 28 | -0.7994486 | 0.4307605 | |
| treatmentHost zero | 8.8224892 | 2.495655 | 28 | 3.5351399 | 0.0014386 |
|
| treatmentMolysis | 4.5783098 | 2.217883 | 28 | 2.0642703 | 0.0483678 |
|
| treatmentQIAamp | 0.8360832 | 2.678401 | 28 | 0.3121575 | 0.7572336 | |
| log10(Final_reads) | 5.6349921 | 1.419898 | 28 | 3.9685903 | 0.0004571 |
|
cat("Species richness at BAL\n\n")
## Species richness at BAL
cat("S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(S.obs ~ treatment + log10 (Final_reads) + (1|original_sample), data = subset(sample_data, sample_data$sample_type == "BAL")) %>%
summary() %>%
.$coefficients %>%
data.frame(check.names = F) %>%
mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Estimate | Std. Error | df | t value | Pr(>|t|) | ||
|---|---|---|---|---|---|---|
| (Intercept) | -100.395549 | 25.504894 | 20.38095 | -3.9363249 | 0.0007925 |
|
| treatmentlyPMA | -5.839094 | 7.201885 | 17.39942 | -0.8107731 | 0.4284473 | |
| treatmentBenzonase | -10.667780 | 7.752673 | 18.42945 | -1.3760131 | 0.1853100 | |
| treatmentHost zero | -8.986746 | 8.091470 | 18.61512 | -1.1106444 | 0.2808621 | |
| treatmentMolysis | 3.342010 | 8.316699 | 18.72132 | 0.4018433 | 0.6923498 | |
| treatmentQIAamp | -10.097992 | 8.350027 | 18.73603 | -1.2093364 | 0.2415732 | |
| log10(Final_reads) | 19.520813 | 4.535759 | 19.97726 | 4.3037587 | 0.0003466 |
|
cat("Species richness at sputum\n\n")
## Species richness at sputum
cat("S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(S.obs ~ treatment + log10 (Final_reads) + (1|original_sample), data = subset(sample_data, sample_data$sample_type == "Sputum")) %>%
summary() %>%
.$coefficients %>%
data.frame(check.names = F) %>%
mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Estimate | Std. Error | df | t value | Pr(>|t|) | ||
|---|---|---|---|---|---|---|
| (Intercept) | -166.73637 | 94.19800 | 20.78908 | -1.770063 | 0.0913826 | |
| treatmentlyPMA | 21.48046 | 13.48708 | 19.68002 | 1.592670 | 0.1271695 | |
| treatmentBenzonase | 41.90046 | 17.06102 | 19.99448 | 2.455918 | 0.0233254 |
|
| treatmentHost zero | 58.57768 | 28.77613 | 20.32431 | 2.035634 | 0.0550337 | |
| treatmentMolysis | 60.65374 | 33.57099 | 20.37080 | 1.806731 | 0.0855997 | |
| treatmentQIAamp | 41.44599 | 24.96239 | 20.26633 | 1.660337 | 0.1122410 | |
| log10(Final_reads) | 31.32813 | 16.05007 | 20.49824 | 1.951899 | 0.0647464 |
#Shannon
cat("Shannon index\n\n")
## Shannon index
cat("Shannon ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## Shannon ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(data_shannon ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample), data = sample_data) %>%
summary() %>%
.$coefficients %>%
data.frame(check.names = F) %>%
mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Estimate | Std. Error | df | t value | Pr(>|t|) | ||
|---|---|---|---|---|---|---|
| (Intercept) | 1.8948018 | 0.8689906 | 73.16148 | 2.1804629 | 0.0324397 |
|
| sample_typeNasal swab | 0.3612177 | 0.4227092 | 52.85894 | 0.8545299 | 0.3966679 | |
| sample_typeSputum | 0.5229775 | 0.4502631 | 46.00549 | 1.1614932 | 0.2514323 | |
| treatmentlyPMA | 0.1432302 | 0.3457153 | 59.27054 | 0.4143009 | 0.6801493 | |
| treatmentBenzonase | 0.3245806 | 0.3510837 | 61.59790 | 0.9245107 | 0.3588282 | |
| treatmentHost zero | 0.1622396 | 0.3591607 | 62.11606 | 0.4517187 | 0.6530450 | |
| treatmentMolysis | 0.9220443 | 0.3646163 | 62.43914 | 2.5288070 | 0.0139853 |
|
| treatmentQIAamp | 0.3416438 | 0.3654292 | 62.48560 | 0.9349109 | 0.3534333 | |
| log10(Final_reads) | -0.2037567 | 0.1481509 | 70.28159 | -1.3753323 | 0.1733965 | |
| sample_typeNasal swab:treatmentlyPMA | -0.3349724 | 0.4621995 | 62.11790 | -0.7247356 | 0.4713354 | |
| sample_typeSputum:treatmentlyPMA | 0.9305639 | 0.4573666 | 58.90852 | 2.0346129 | 0.0463989 |
|
| sample_typeNasal swab:treatmentBenzonase | -0.1716183 | 0.4431174 | 62.20574 | -0.3872975 | 0.6998579 | |
| sample_typeSputum:treatmentBenzonase | 1.1940106 | 0.4466493 | 59.19803 | 2.6732618 | 0.0096921 |
|
| sample_typeNasal swab:treatmentHost zero | 0.1836222 | 0.4309190 | 61.71186 | 0.4261176 | 0.6715048 | |
| sample_typeSputum:treatmentHost zero | 1.4844081 | 0.4573881 | 59.56239 | 3.2454017 | 0.0019263 |
|
| sample_typeNasal swab:treatmentMolysis | -0.6161157 | 0.4488821 | 62.42148 | -1.3725557 | 0.1748022 | |
| sample_typeSputum:treatmentMolysis | 0.9186726 | 0.4659713 | 59.93997 | 1.9715216 | 0.0532853 | |
| sample_typeNasal swab:treatmentQIAamp | -0.2466370 | 0.4305635 | 61.76297 | -0.5728238 | 0.5688436 | |
| sample_typeSputum:treatmentQIAamp | 1.0739999 | 0.4489918 | 59.21646 | 2.3920256 | 0.0199485 |
|
#Simpson
cat("Inverse Simpson index\n\n")
## Inverse Simpson index
cat("InvSimp ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## InvSimp ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(data_invsimpson ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample), data = sample_data) %>%
summary() %>%
.$coefficients %>%
data.frame(check.names = F) %>%
mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Estimate | Std. Error | df | t value | Pr(>|t|) | ||
|---|---|---|---|---|---|---|
| (Intercept) | 5.0016985 | 4.1346231 | 73.75264 | 1.2097109 | 0.2302531 | |
| sample_typeNasal swab | -0.1347032 | 1.9113437 | 60.99582 | -0.0704757 | 0.9440456 | |
| sample_typeSputum | 0.3376540 | 2.0144643 | 55.69325 | 0.1676148 | 0.8674935 | |
| treatmentlyPMA | -0.0446283 | 1.6994301 | 60.26753 | -0.0262607 | 0.9791362 | |
| treatmentBenzonase | -0.0159936 | 1.7179035 | 63.20493 | -0.0093099 | 0.9926012 | |
| treatmentHost zero | -0.5969192 | 1.7555995 | 63.84335 | -0.3400088 | 0.7349661 | |
| treatmentMolysis | 2.8388719 | 1.7811019 | 64.24055 | 1.5938852 | 0.1158705 | |
| treatmentQIAamp | -0.2830161 | 1.7849044 | 64.29760 | -0.1585609 | 0.8745118 | |
| log10(Final_reads) | -0.4205707 | 0.7102118 | 72.89587 | -0.5921764 | 0.5555648 | |
| sample_typeNasal swab:treatmentlyPMA | -0.2501880 | 2.2597896 | 63.56198 | -0.1107130 | 0.9121928 | |
| sample_typeSputum:treatmentlyPMA | 3.6784976 | 2.2498385 | 59.81795 | 1.6350052 | 0.1073001 | |
| sample_typeNasal swab:treatmentBenzonase | 0.5135753 | 2.1660273 | 63.71628 | 0.2371047 | 0.8133362 | |
| sample_typeSputum:treatmentBenzonase | 8.4907847 | 2.1958416 | 60.21080 | 3.8667565 | 0.0002731 |
|
| sample_typeNasal swab:treatmentHost zero | 1.1519699 | 2.1086029 | 63.04119 | 0.5463190 | 0.5867750 | |
| sample_typeSputum:treatmentHost zero | 7.4853477 | 2.2470527 | 60.66550 | 3.3311847 | 0.0014760 |
|
| sample_typeNasal swab:treatmentMolysis | -2.0541758 | 2.1931699 | 64.02401 | -0.9366241 | 0.3524716 | |
| sample_typeSputum:treatmentMolysis | 4.7103881 | 2.2875456 | 61.13466 | 2.0591450 | 0.0437470 |
|
| sample_typeNasal swab:treatmentQIAamp | 0.4812870 | 2.1067028 | 63.07313 | 0.2284551 | 0.8200315 | |
| sample_typeSputum:treatmentQIAamp | 5.7791599 | 2.2072771 | 60.23515 | 2.6182303 | 0.0111632 |
|
#BPI
cat("Beger Parker index\n\n")
## Beger Parker index
cat("BPI ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## BPI ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(dbp ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample), data = sample_data) %>%
summary() %>%
.$coefficients %>%
data.frame(check.names = F) %>%
mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Estimate | Std. Error | df | t value | Pr(>|t|) | ||
|---|---|---|---|---|---|---|
| (Intercept) | -0.0429030 | 0.3089340 | 73.34596 | -0.1388743 | 0.8899301 | |
| sample_typeNasal swab | -0.2254845 | 0.1470717 | 55.46609 | -1.5331600 | 0.1309210 | |
| sample_typeSputum | -0.1615566 | 0.1559725 | 48.85446 | -1.0358016 | 0.3053947 | |
| treatmentlyPMA | -0.0992353 | 0.1244133 | 59.03815 | -0.7976258 | 0.4282842 | |
| treatmentBenzonase | -0.2090856 | 0.1261407 | 61.68022 | -1.6575576 | 0.1024850 | |
| treatmentHost zero | -0.1071198 | 0.1289953 | 62.26501 | -0.8304161 | 0.4094764 | |
| treatmentMolysis | -0.3789439 | 0.1309245 | 62.62965 | -2.8943694 | 0.0052227 |
|
| treatmentQIAamp | -0.2002048 | 0.1312121 | 62.68208 | -1.5258106 | 0.1320873 | |
| log10(Final_reads) | 0.1455197 | 0.0528497 | 71.25889 | 2.7534627 | 0.0074764 |
|
| sample_typeNasal swab:treatmentlyPMA | 0.1074786 | 0.1660074 | 62.17039 | 0.6474324 | 0.5197358 | |
| sample_typeSputum:treatmentlyPMA | -0.1999551 | 0.1646342 | 58.63042 | -1.2145416 | 0.2294101 | |
| sample_typeNasal swab:treatmentBenzonase | 0.0783117 | 0.1591430 | 62.28479 | 0.4920838 | 0.6243901 | |
| sample_typeSputum:treatmentBenzonase | -0.2242934 | 0.1607441 | 58.96684 | -1.3953444 | 0.1681450 | |
| sample_typeNasal swab:treatmentHost zero | -0.0843810 | 0.1548174 | 61.70491 | -0.5450359 | 0.5876935 | |
| sample_typeSputum:treatmentHost zero | -0.4349677 | 0.1645678 | 59.37801 | -2.6430911 | 0.0104899 |
|
| sample_typeNasal swab:treatmentMolysis | 0.1755526 | 0.1611877 | 62.54326 | 1.0891189 | 0.2802803 | |
| sample_typeSputum:treatmentMolysis | -0.2575968 | 0.1676125 | 59.80380 | -1.5368593 | 0.1296021 | |
| sample_typeNasal swab:treatmentQIAamp | 0.1096753 | 0.1546846 | 61.75206 | 0.7090255 | 0.4809770 | |
| sample_typeSputum:treatmentQIAamp | -0.2873906 | 0.1615851 | 58.98802 | -1.7785712 | 0.0804633 |
A4. Taxa beta diversity
Permanova (Taxa dist ~ log10(final reads) + sample_type + treatment + sample_type * treatment + subject_id) –> both stratified and nonstratified
phyloseq_rel_nz <- subset_samples(phyloseq$phyloseq_rel, S.obs != 0 & sample_type %in% c("BAL", "Nasal swab", "Sputum"))
bray_perm_uni <- vegan::adonis2(distance(phyloseq_rel_nz, method="bray") ~ sample_type + log10(Final_reads) + treatment + subject_id,
data = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F), permutations = 10000)
bray_perm <- vegan::adonis2(distance(phyloseq_rel_nz, method="bray") ~ sample_type + log10(Final_reads) + lypma + benzonase + host_zero + molysis + qiaamp + subject_id,
data = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F), permutations = 10000)
bray_perm_strata <- vegan::adonis2(distance(phyloseq_rel_nz, method="bray") ~ sample_type + log10(Final_reads) + lypma + benzonase + host_zero + molysis + qiaamp,
data = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F),
strata = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F) %>% .$subject_id, permutations = 10000)
bray_perm_inter <- vegan::adonis2(distance(phyloseq_rel_nz, method="bray") ~ sample_type * treatment + log10(Final_reads) + subject_id,
data = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F), permutations = 10000)
bray_perm_ns <- vegan::adonis2(distance(subset_samples(phyloseq_rel_nz, sample_type == "Nasal swab"), method="bray") ~ lypma + benzonase + host_zero + molysis + qiaamp + log10(Final_reads) + subject_id,
data = subset_samples(phyloseq_rel_nz, sample_type == "Nasal swab") %>% sample_data %>% data.frame(check.names = F), permutations = 10000)
bray_perm_bal <- vegan::adonis2(distance(subset_samples(phyloseq_rel_nz, sample_type == "BAL"), method="bray") ~ lypma + benzonase + host_zero + molysis + qiaamp + log10(Final_reads) + subject_id,
data = subset_samples(phyloseq_rel_nz, sample_type == "BAL") %>% sample_data %>% data.frame(check.names = F), permutations = 10000)
bray_perm_spt <- vegan::adonis2(distance(subset_samples(phyloseq_rel_nz, sample_type == "Sputum"), method="bray") ~ lypma + benzonase + host_zero + molysis + qiaamp + log10(Final_reads) + subject_id,
data = subset_samples(phyloseq_rel_nz, sample_type == "Sputum") %>% sample_data %>% data.frame(check.names = F), permutations = 10000)
cat("\nUnivariate analysis\n")
##
## Univariate analysis
bray_perm_uni %>% data.frame(check.names = F) %>% rownames_to_column("row.names") %>%
mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
row.names == "treatment" ~ 'Treatment',
row.names == "subject_id" ~ 'Subject',
row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
row.names == "Residual" ~ 'Residual',
row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>%
round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Df | SumOfSqs | R2 | F | Pr(>F) | ||
|---|---|---|---|---|---|---|
| Sample type | 2 | 13.429 | 0.371 | 56.251 | 0.000 |
|
| log10(Final reads) | 1 | 1.111 | 0.031 | 9.310 | 0.000 |
|
| Treatment | 5 | 1.213 | 0.033 | 2.033 | 0.001 |
|
| Subject | 10 | 11.639 | 0.321 | 9.751 | 0.000 |
|
| Residual | 74 | 8.833 | 0.244 | NA | NA | |
| Total | 92 | 36.225 | 1.000 | NA | NA |
cat("\nDetailed treatment\n")
##
## Detailed treatment
bray_perm %>% data.frame(check.names = F) %>% rownames_to_column("row.names") %>%
mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
row.names == "lypma" ~ 'lyPMA',
row.names == "benzonase" ~ 'Benzonase',
row.names == "host_zero" ~ 'Host zero',
row.names == "molysis" ~ 'Molysis',
row.names == "qiaamp" ~ 'QIAamp',
row.names == "subject_id" ~ 'Subject',
row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
row.names == "Residual" ~ 'Residual',
row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>%
round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Df | SumOfSqs | R2 | F | Pr(>F) | ||
|---|---|---|---|---|---|---|
| Sample type | 2 | 13.429 | 0.371 | 56.251 | 0.000 |
|
| log10(Final reads) | 1 | 1.111 | 0.031 | 9.310 | 0.000 |
|
| lyPMA | 1 | 0.140 | 0.004 | 1.172 | 0.287 | |
| Benzonase | 1 | 0.092 | 0.003 | 0.774 | 0.607 | |
| Host zero | 1 | 0.157 | 0.004 | 1.312 | 0.218 | |
| Molysis | 1 | 0.200 | 0.006 | 1.676 | 0.103 | |
| QIAamp | 1 | 0.624 | 0.017 | 5.229 | 0.000 |
|
| Subject | 10 | 11.639 | 0.321 | 9.751 | 0.000 |
|
| Residual | 74 | 8.833 | 0.244 | NA | NA | |
| Total | 92 | 36.225 | 1.000 | NA | NA |
cat("\n Strata -detailed treatment\n")
##
## Strata -detailed treatment
bray_perm_strata %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>%
mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
row.names == "lypma" ~ 'lyPMA',
row.names == "benzonase" ~ 'Benzonase',
row.names == "host_zero" ~ 'Host zero',
row.names == "molysis" ~ 'Molysis',
row.names == "qiaamp" ~ 'QIAamp',
row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
row.names == "Residual" ~ 'Residual',
row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>%
round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Df | SumOfSqs | R2 | F | Pr(>F) | ||
|---|---|---|---|---|---|---|
| Sample type | 2 | 13.429 | 0.371 | 27.550 | 0.000 |
|
| log10(Final reads) | 1 | 1.111 | 0.031 | 4.560 | 0.000 |
|
| lyPMA | 1 | 0.140 | 0.004 | 0.574 | 0.530 | |
| Benzonase | 1 | 0.092 | 0.003 | 0.379 | 0.665 | |
| Host zero | 1 | 0.157 | 0.004 | 0.643 | 0.406 | |
| Molysis | 1 | 0.200 | 0.006 | 0.821 | 0.167 | |
| QIAamp | 1 | 0.624 | 0.017 | 2.561 | 0.003 |
|
| Residual | 84 | 20.472 | 0.565 | NA | NA | |
| Total | 92 | 36.225 | 1.000 | NA | NA |
cat("\n Interaction term \n")
##
## Interaction term
bray_perm_inter %>% data.frame(check.names = F) %>% rownames_to_column("row.names") %>%
mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
row.names == "treatment" ~ 'Treatment',
row.names == "subject_id" ~ 'Subject',
row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
row.names == "sample_type:treatment" ~ 'Sample type X treatment',
row.names == "Residual" ~ 'Residual',
row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>%
round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Df | SumOfSqs | R2 | F | Pr(>F) | ||
|---|---|---|---|---|---|---|
| Sample type | 2 | 13.429 | 0.371 | 65.605 | 0 |
|
| Treatment | 5 | 1.184 | 0.033 | 2.314 | 0 |
|
| log10(Final reads) | 1 | 1.140 | 0.031 | 11.142 | 0 |
|
| Subject | 10 | 11.639 | 0.321 | 11.372 | 0 |
|
| Sample type X treatment | 10 | 2.283 | 0.063 | 2.231 | 0 |
|
| Residual | 64 | 6.550 | 0.181 | NA | NA | |
| Total | 92 | 36.225 | 1.000 | NA | NA |
cat("\n Stratified - nasal swab \n")
##
## Stratified - nasal swab
bray_perm_ns %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>%
mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
row.names == "benzonase" ~ 'Benzonase',
row.names == "host_zero" ~ 'Host zero',
row.names == "molysis" ~ 'Molysis',
row.names == "qiaamp" ~ 'QIAamp',
row.names == "subject_id" ~ 'Subject id',
row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
row.names == "Residual" ~ 'Residual',
row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>%
round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Df | SumOfSqs | R2 | F | Pr(>F) | ||
|---|---|---|---|---|---|---|
| lyPMA | 1 | 0.730 | 0.120 | 5.122 | 0.004 |
|
| Benzonase | 1 | 0.191 | 0.031 | 1.342 | 0.262 | |
| Host zero | 1 | 0.171 | 0.028 | 1.201 | 0.306 | |
| Molysis | 1 | 0.137 | 0.022 | 0.961 | 0.405 | |
| QIAamp | 1 | 0.254 | 0.042 | 1.780 | 0.148 | |
| log10(Final reads) | 1 | 0.428 | 0.070 | 3.007 | 0.036 |
|
| Subject id | 2 | 0.488 | 0.080 | 1.714 | 0.119 | |
| Residual | 26 | 3.704 | 0.607 | NA | NA | |
| Total | 34 | 6.103 | 1.000 | NA | NA |
cat("\n Stratified - BAL \n")
##
## Stratified - BAL
bray_perm_bal %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>%
mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
row.names == "benzonase" ~ 'Benzonase',
row.names == "host_zero" ~ 'Host zero',
row.names == "molysis" ~ 'Molysis',
row.names == "qiaamp" ~ 'QIAamp',
row.names == "subject_id" ~ 'Subject id',
row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
row.names == "Residual" ~ 'Residual',
row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>%
round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Df | SumOfSqs | R2 | F | Pr(>F) | ||
|---|---|---|---|---|---|---|
| lyPMA | 1 | 0.100 | 0.010 | 1.144 | 0.308 | |
| Benzonase | 1 | 0.025 | 0.003 | 0.288 | 0.934 | |
| Host zero | 1 | 0.086 | 0.009 | 0.990 | 0.418 | |
| Molysis | 1 | 0.085 | 0.009 | 0.969 | 0.428 | |
| QIAamp | 1 | 0.229 | 0.024 | 2.625 | 0.035 |
|
| log10(Final reads) | 1 | 1.482 | 0.152 | 16.966 | 0.000 |
|
| Subject id | 4 | 6.241 | 0.641 | 17.864 | 0.000 |
|
| Residual | 17 | 1.485 | 0.153 | NA | NA | |
| Total | 27 | 9.734 | 1.000 | NA | NA |
cat("\n Stratified - sputum \n")
##
## Stratified - sputum
bray_perm_spt %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>%
mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
row.names == "benzonase" ~ 'Benzonase',
row.names == "host_zero" ~ 'Host zero',
row.names == "molysis" ~ 'Molysis',
row.names == "qiaamp" ~ 'QIAamp',
row.names == "subject_id" ~ 'Subject id',
row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
row.names == "Residual" ~ 'Residual',
row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>%
round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Df | SumOfSqs | R2 | F | Pr(>F) | ||
|---|---|---|---|---|---|---|
| lyPMA | 1 | 0.139 | 0.020 | 2.565 | 0.046 |
|
| Benzonase | 1 | 0.037 | 0.005 | 0.688 | 0.635 | |
| Host zero | 1 | 0.171 | 0.025 | 3.150 | 0.018 |
|
| Molysis | 1 | 0.436 | 0.063 | 8.046 | 0.000 |
|
| QIAamp | 1 | 0.953 | 0.137 | 17.588 | 0.000 |
|
| log10(Final reads) | 1 | 0.172 | 0.025 | 3.174 | 0.017 |
|
| Subject id | 4 | 4.022 | 0.578 | 18.555 | 0.000 |
|
| Residual | 19 | 1.030 | 0.148 | NA | NA | |
| Total | 29 | 6.960 | 1.000 | NA | NA |
A5. DA analysis for taxa, by sample type and treatment
--> both stratified and nonstratified
--> look at level groups as well - phylum, etc.
<<<<<<< HEAD
#DA analysis - MaAslin
sample_data(phyloseq_rel_nz)$log10.Final_reads <- log10(sample_data(phyloseq_rel_nz)$Final_reads)
sample_data(phyloseq_rel_nz)$sampletype_treatment <- paste(sample_data(phyloseq_rel_nz)$sample_type, sample_data(phyloseq_rel_nz)$treatment, sep = ":")
#Running MaAslin for all sample without decontam
#for taxa differentially abundant by host depletion method, look to see which ones overlap with potential contaminant taxa
# Maaslin - # # y ~ log(final reads) + sample_type + treatment -----------
#all samples
maaslin_all = Maaslin2(input_data = otu_table(phyloseq_rel_nz) %>% t %>% data.frame(),
input_metadata = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F),
output = "/Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output",
fixed_effects = c("sample_type", "log10.Final_reads","lypma", "benzonase", "host_zero", "molysis", "qiaamp"),
transform = "LOG", #default
normalization = "TSS", # default
random_effects = c("subject_id"),
reference = c("sample_type,BAL"),
plot_heatmap = F,
plot_scatter = F)
## [1] "Warning: Deleting existing log file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/maaslin2.log"
## 2023-03-22 00:45:26 INFO::Writing function arguments to log file
## 2023-03-22 00:45:26 INFO::Verifying options selected are valid
## 2023-03-22 00:45:26 INFO::Determining format of input files
## 2023-03-22 00:45:26 INFO::Input format is data samples as rows and metadata samples as rows
## 2023-03-22 00:45:26 INFO::Formula for random effects: expr ~ (1 | subject_id)
## 2023-03-22 00:45:26 INFO::Formula for fixed effects: expr ~ sample_type + log10.Final_reads + lypma + benzonase + host_zero + molysis + qiaamp
## 2023-03-22 00:45:26 INFO::Filter data based on min abundance and min prevalence
## 2023-03-22 00:45:26 INFO::Total samples in data: 93
## 2023-03-22 00:45:26 INFO::Min samples required with min abundance for a feature not to be filtered: 9.300000
## 2023-03-22 00:45:26 INFO::Total filtered features: 183
## 2023-03-22 00:45:26 INFO::Filtered feature names from abundance and prevalence filtering: Actinomyces_sp_oral_taxon_414, Actinomyces_sp_oral_taxon_448, Actinomyces_sp_oral_taxon_897, Parascardovia_denticolens, Corynebacterium_afermentans, Corynebacterium_aurimucosum, Corynebacterium_kroppenstedtii, Corynebacterium_matruchotii, Corynebacterium_tuberculostearicum, Mycobacterium_chimaera, Mycobacterium_intracellulare, Micrococcus_luteus, Tropheryma_whipplei, Cutibacterium_avidum, Atopobium_deltae, Adlercreutzia_equolifaciens, Rubrobacter_radiotolerans, Patulibacter_medicamentivorans, Prevotella_oris, Capnocytophaga_leadbetteri, Brochothrix_campestris, Staphylococcus_capitis, Enterococcus_faecalis, Lactobacillus_rhamnosus, Streptococcus_sinensis, Eubacterium_nodatum, Lachnospiraceae_bacterium_oral_taxon_096, Oribacterium_parvum, Shuttleworthia_satelles, Selenomonas_noxia, Selenomonas_sputigena, Anaerococcus_nagyae, Anaerococcus_octavius, Peptoniphilus_lacrimalis, Leptotrichia_wadei, Achromobacter_denitrificans, Achromobacter_ruhlandii, Achromobacter_xylosoxidans, Lautropia_mirabilis, Leptothrix_ochracea, Cardiobacterium_hominis, Escherichia_coli, Klebsiella_michiganensis, Klebsiella_oxytoca, Serratia_liquefaciens, Fretibacterium_fastidiosum, Aspergillus_eucalypticola, Aspergillus_kawachii, Aspergillus_lacticoffeatus, Aspergillus_niger, Aspergillus_phoenicis, Aspergillus_sydowii, Aspergillus_thermomutatus, Aspergillus_tubingensis, Aspergillus_turcosus, Aspergillus_vadensis, Aspergillus_welwitschiae, Candida_tropicalis, Malassezia_globosa, Actinomyces_cardiffensis, Actinomyces_denticolens, Actinomyces_radingae, Actinomyces_turicensis, Varibaculum_cambriense, Bifidobacterium_breve, Gardnerella_vaginalis, Scardovia_inopinata, Corynebacterium_pyruviciproducens, Gordonia_polyisoprenivorans, Mycolicibacterium_fortuitum, Mycolicibacterium_rufum, Nakamurella_silvestris, Atopobium_minutum, Gordonibacter_pamelaeae, Bacteroidetes_oral_taxon_274, Porphyromonas_asaccharolytica, Porphyromonas_canoris, Porphyromonas_catoniae, Porphyromonas_uenonis, Alloprevotella_rava, Alloprevotella_tannerae, Prevotella_buccae, Prevotella_denticola, Prevotella_intermedia, Prevotella_nigrescens, Prevotella_oulorum, Prevotella_pleuritidis, Prevotella_scopos, Prevotella_sp_F0091, Prevotella_sp_oral_taxon_306, Prevotella_veroralis, Tannerella_forsythia, Tannerella_sp_oral_taxon_808, Capnocytophaga_granulosa, Capnocytophaga_sputigena, Litorilinea_aerophila, Bacillus_cereus_group, Bacillus_ginsengihumi, Listeria_floridensis, Agitococcus_lubricus, Lactobacillus_fermentum, Lactobacillus_gasseri, Lactobacillus_oris, Lactobacillus_paragasseri, Lactobacillus_reuteri, Lactobacillus_salivarius, Sharpea_azabuensis, Streptococcus_massiliensis, Streptococcus_sp_SK643, Streptococcus_sp_oral_taxon_056, Mageeibacillus_indolicus, Oribacterium_asaccharolyticum, Peptostreptococcus_sp_MV1, Leptotrichia_sp_oral_taxon_212, Desulfobulbus_oralis, Moraxella_catarrhalis, Stenotrophomonas_pavanii, Treponema_socranskii, Bacillus_intestinalis, Listeria_innocua, Listeria_monocytogenes, Staphylococcus_haemolyticus, Enterococcus_avium, Streptococcus_downei, Streptococcus_pyogenes, Streptococcus_salivarius_CAG_79, Streptococcus_sobrinus, Streptococcus_sp_DD11, Streptococcus_sp_HMSC070B10, Streptococcus_sp_NLAE_zl_C503, Streptococcus_thermophilus, Streptococcus_viridans, Eubacterium_saphenum, Eubacterium_rectale, Filifactor_alocis, Eggerthia_catenaformis, Selenomonas_flueggei, Selenomonas_sp_oral_taxon_892, Selenomonas_sp_oral_taxon_920, Dialister_invisus, Dialister_micraerophilus, Negativicoccus_succinicivorans, Veillonella_tobetsuensis, Fusobacterium_periodonticum, Fusobacterium_sp_oral_taxon_370, Leptotrichia_buccalis, Leptotrichia_hofstadii, Leptotrichia_sp_oral_taxon_215, Leptotrichia_sp_oral_taxon_225, Leptotrichia_sp_oral_taxon_498, Leptotrichia_sp_oral_taxon_879, Fimbriiglobus_ruber, Rickettsia_felis, Burkholderia_multivorans, Eikenella_corrodens, Neisseria_elongata, Campylobacter_showae, Haemophilus_parainfluenzae, Candida_orthopsilosis, Rickettsia_typhi, Neisseria_macacae, Neisseria_mucosa, Neisseria_sicca, Nannocystis_exedens, Campylobacter_concisus, Campylobacter_gracilis, Campylobacter_mucosalis, Franconibacter_helveticus, Klebsiella_pneumoniae, Salmonella_enterica, Superficieibacter_electus, Serratia_marcescens, Haemophilus_sp_HMSC71H05, Pseudomonas_fluorescens_group, Pseudomonas_geniculata, Aspergillus_fumigatus, Saccharomyces_cerevisiae, Saccharomyces_cerevisiae_x_Saccharomyces_kudriavzevii, Saccharomyces_kudriavzevii, Cryptococcus_gattii_VGI, Cryptococcus_gattii_VGII, Cryptococcus_gattii_VGIII, Cryptococcus_neoformans
## 2023-03-22 00:45:26 INFO::Total filtered features with variance filtering: 0
## 2023-03-22 00:45:26 INFO::Filtered feature names from variance filtering:
## 2023-03-22 00:45:26 INFO::Running selected normalization method: TSS
## 2023-03-22 00:45:26 INFO::Applying z-score to standardize continuous metadata
## 2023-03-22 00:45:26 INFO::Running selected transform method: LOG
## 2023-03-22 00:45:26 INFO::Running selected analysis method: LM
## 2023-03-22 00:45:26 INFO::Fitting model to feature number 1, Actinobaculum_sp_oral_taxon_183
## 2023-03-22 00:45:26 INFO::Fitting model to feature number 2, Actinomyces_georgiae
## 2023-03-22 00:45:26 INFO::Fitting model to feature number 3, Actinomyces_graevenitzii
## 2023-03-22 00:45:26 INFO::Fitting model to feature number 4, Actinomyces_hongkongensis
## 2023-03-22 00:45:26 INFO::Fitting model to feature number 5, Actinomyces_johnsonii
## 2023-03-22 00:45:26 INFO::Fitting model to feature number 6, Actinomyces_massiliensis
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 7, Actinomyces_meyeri
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 8, Actinomyces_naeslundii
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 9, Actinomyces_odontolyticus
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 10, Actinomyces_oris
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 11, Actinomyces_sp_HMSC035G02
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 12, Actinomyces_sp_HPA0247
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 13, Actinomyces_sp_ICM47
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 14, Actinomyces_sp_S6_Spd3
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 15, Actinomyces_sp_oral_taxon_170
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 16, Actinomyces_sp_oral_taxon_180
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 17, Actinomyces_sp_oral_taxon_181
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 18, Actinomyces_viscosus
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 19, Aeriscardovia_aeriphila
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 20, Alloscardovia_omnicolens
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 21, Bifidobacterium_dentium
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 22, Bifidobacterium_longum
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 23, Scardovia_wiggsiae
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 24, Corynebacterium_accolens
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 25, Corynebacterium_atypicum
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 26, Corynebacterium_durum
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 27, Corynebacterium_pseudodiphtheriticum
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 28, Corynebacterium_pseudogenitalium
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 29, Rothia_aeria
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 30, Rothia_dentocariosa
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 31, Rothia_mucilaginosa
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 32, Cutibacterium_acnes
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 33, Cutibacterium_granulosum
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 34, Propionibacterium_humerusii
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 35, Propionibacterium_namnetense
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 36, Propionibacterium_acidifaciens
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 37, Pseudopropionibacterium_propionicum
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 38, Atopobium_parvulum
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 39, Atopobium_rimae
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 40, Olsenella_profusa
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 41, Olsenella_scatoligenes
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 42, Olsenella_uli
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 43, Collinsella_aerofaciens
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 44, Collinsella_intestinalis
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 45, Collinsella_stercoris
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 46, Enorma_massiliensis
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 47, X.Collinsella._massiliensis
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 48, Asaccharobacter_celatus
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 49, Denitrobacterium_detoxificans
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 50, Enterorhabdus_caecimuris
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 51, Slackia_exigua
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 52, Slackia_isoflavoniconvertens
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 53, Thermoleophilum_album
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 54, Porphyromonas_endodontalis
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 55, Porphyromonas_somerae
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 56, Prevotella_histicola
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 57, Prevotella_melaninogenica
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 58, Kouleothrix_aurantiaca
## 2023-03-22 00:45:27 INFO::Fitting model to feature number 59, Hydrogenibacillus_schlegelii
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 60, Gemella_asaccharolytica
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 61, Gemella_bergeri
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 62, Gemella_haemolysans
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 63, Gemella_morbillorum
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 64, Gemella_sanguinis
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 65, Staphylococcus_argenteus
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 66, Staphylococcus_aureus
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 67, Staphylococcus_epidermidis
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 68, Staphylococcus_schweitzeri
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 69, Abiotrophia_defectiva
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 70, Abiotrophia_sp_HMSC24B09
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 71, Dolosigranulum_pigrum
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 72, Granulicatella_adiacens
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 73, Granulicatella_elegans
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 74, Streptococcus_anginosus_group
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 75, Streptococcus_australis
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 76, Streptococcus_cristatus
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 77, Streptococcus_gordonii
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 78, Streptococcus_infantis
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 79, Streptococcus_milleri
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 80, Streptococcus_mitis
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 81, Streptococcus_oralis
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 82, Streptococcus_parasanguinis
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 83, Streptococcus_peroris
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 84, Streptococcus_pneumoniae
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 85, Streptococcus_pseudopneumoniae
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 86, Streptococcus_salivarius
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 87, Streptococcus_sanguinis
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 88, Streptococcus_sp_A12
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 89, Streptococcus_sp_F0442
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 90, Streptococcus_sp_HMSC034E03
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 91, Streptococcus_sp_HMSC067H01
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 92, Streptococcus_sp_HMSC071D03
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 93, Streptococcus_sp_HPH0090
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 94, Streptococcus_sp_M334
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 95, Streptococcus_vestibularis
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 96, Eubacterium_brachy
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 97, Eubacterium_infirmum
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 98, Eubacterium_sulci
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 99, Mogibacterium_diversum
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 100, Mogibacterium_pumilum
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 101, Mogibacterium_timidum
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 102, Lachnoanaerobaculum_saburreum
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 103, Oribacterium_sinus
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 104, Oribacterium_sp_oral_taxon_078
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 105, Stomatobaculum_longum
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 106, Peptostreptococcus_stomatis
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 107, Bulleidia_extructa
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 108, Solobacterium_moorei
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 109, Limnochorda_pilosa
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 110, Veillonella_atypica
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 111, Veillonella_dispar
## 2023-03-22 00:45:28 INFO::Fitting model to feature number 112, Veillonella_infantium
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 113, Veillonella_parvula
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 114, Veillonella_sp_T11011_6
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 115, Finegoldia_magna
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 116, Parvimonas_micra
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 117, Parvimonas_sp_oral_taxon_110
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 118, Parvimonas_sp_oral_taxon_393
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 119, Fusobacterium_nucleatum
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 120, Gemmata_obscuriglobus
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 121, Paludisphaera_borealis
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 122, Cupriavidus_sp
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 123, Sutterella_parvirubra
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 124, Anaerobiospirillum_thomasii
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 125, Cardiobacterium_valvarum
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 126, Alkalilimnicola_ehrlichii
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 127, Thiohalorhabdus_denitrificans
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 128, Pseudomonas_aeruginosa_group
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 129, Acholeplasma_oculi
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 130, Candida_albicans
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 131, Candida_dubliniensis
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 132, Candida_parapsilosis
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 133, Malassezia_restricta
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 134, Prevotella_jejuni
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 135, Prevotella_pallens
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 136, Prevotella_salivae
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 137, Tannerella_sp_oral_taxon_HOT_286
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 138, Capnocytophaga_gingivalis
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 139, Streptococcus_mutans
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 140, Megasphaera_micronuciformis
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 141, Neisseria_flavescens
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 142, Neisseria_subflava
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 143, Stenotrophomonas_maltophilia
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 144, Stenotrophomonas_rhizophila
## 2023-03-22 00:45:29 INFO::Counting total values for each feature
## 2023-03-22 00:45:29 WARNING::Deleting existing residuals file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/residuals.rds
## 2023-03-22 00:45:29 INFO::Writing residuals to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/residuals.rds
## 2023-03-22 00:45:29 WARNING::Deleting existing fitted file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/fitted.rds
## 2023-03-22 00:45:29 INFO::Writing fitted values to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/fitted.rds
## 2023-03-22 00:45:29 WARNING::Deleting existing ranef file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/ranef.rds
## 2023-03-22 00:45:29 INFO::Writing extracted random effects to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/ranef.rds
## 2023-03-22 00:45:29 INFO::Writing all results to file (ordered by increasing q-values): /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/all_results.tsv
## 2023-03-22 00:45:29 INFO::Writing the significant results (those which are less than or equal to the threshold of 0.250000 ) to file (ordered by increasing q-values): /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/significant_results.tsv
#Checking number of bugs for the main text of the manuscript
cat("Number of differentially abundant bugs by each metadata")
## Number of differentially abundant bugs by each metadata
maaslin_all$results %>% subset(., .$qval < 0.1) %>% .$metadata %>% table()
## .
## benzonase host_zero log10.Final_reads lypma
## 8 6 50 14
## molysis qiaamp sample_type
## 9 2 66
cat("Decreased bugs by each metadata")
## Decreased bugs by each metadata
maaslin_all$results %>% subset(., .$qval < 0.1 & .$coef < 0) %>% .$metadata %>% table()
## .
## benzonase host_zero log10.Final_reads lypma
## 1 1 1 2
## molysis sample_type
## 1 1
cat("Incrased")
## Incrased
maaslin_all$results %>% subset(., .$qval < 0.1 & .$coef > 0) %>% .$metadata %>% table()
## .
## benzonase host_zero log10.Final_reads lypma
## 7 5 49 12
## molysis qiaamp sample_type
## 8 2 65
#NS
# # y ~ log(final reads) + sample_type + treatment
fit_data_ns = Maaslin2(input_data = otu_table(subset_samples(phyloseq_rel_nz, sample_type == "Nasal swab")) %>% t %>% data.frame(),
input_metadata = sample_data(subset_samples(phyloseq_rel_nz, sample_type == "Nasal swab")) %>% data.frame(),
output = "/Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output",
fixed_effects = c("log10.Final_reads","lypma", "benzonase", "host_zero", "molysis", "qiaamp"),
transform = "LOG", #default
normalization = "TSS", # default
random_effects = c("subject_id"),
plot_heatmap = F,
plot_scatter = F)
## [1] "Warning: Deleting existing log file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/maaslin2.log"
## 2023-03-22 00:45:29 INFO::Writing function arguments to log file
## 2023-03-22 00:45:29 INFO::Verifying options selected are valid
## 2023-03-22 00:45:29 INFO::Determining format of input files
## 2023-03-22 00:45:29 INFO::Input format is data samples as rows and metadata samples as rows
## 2023-03-22 00:45:29 INFO::Formula for random effects: expr ~ (1 | subject_id)
## 2023-03-22 00:45:29 INFO::Formula for fixed effects: expr ~ log10.Final_reads + lypma + benzonase + host_zero + molysis + qiaamp
## 2023-03-22 00:45:29 INFO::Filter data based on min abundance and min prevalence
## 2023-03-22 00:45:29 INFO::Total samples in data: 35
## 2023-03-22 00:45:29 INFO::Min samples required with min abundance for a feature not to be filtered: 3.500000
## 2023-03-22 00:45:29 INFO::Total filtered features: 295
## 2023-03-22 00:45:29 INFO::Filtered feature names from abundance and prevalence filtering: Actinobaculum_sp_oral_taxon_183, Actinomyces_georgiae, Actinomyces_graevenitzii, Actinomyces_hongkongensis, Actinomyces_johnsonii, Actinomyces_massiliensis, Actinomyces_meyeri, Actinomyces_naeslundii, Actinomyces_odontolyticus, Actinomyces_oris, Actinomyces_sp_HMSC035G02, Actinomyces_sp_HPA0247, Actinomyces_sp_ICM47, Actinomyces_sp_S6_Spd3, Actinomyces_sp_oral_taxon_170, Actinomyces_sp_oral_taxon_180, Actinomyces_sp_oral_taxon_181, Actinomyces_sp_oral_taxon_414, Actinomyces_sp_oral_taxon_448, Actinomyces_sp_oral_taxon_897, Actinomyces_viscosus, Alloscardovia_omnicolens, Bifidobacterium_dentium, Bifidobacterium_longum, Parascardovia_denticolens, Scardovia_wiggsiae, Corynebacterium_afermentans, Corynebacterium_durum, Corynebacterium_kroppenstedtii, Corynebacterium_matruchotii, Corynebacterium_tuberculostearicum, Mycobacterium_chimaera, Mycobacterium_intracellulare, Micrococcus_luteus, Rothia_aeria, Rothia_dentocariosa, Tropheryma_whipplei, Propionibacterium_acidifaciens, Pseudopropionibacterium_propionicum, Atopobium_deltae, Atopobium_parvulum, Atopobium_rimae, Olsenella_profusa, Olsenella_scatoligenes, Olsenella_uli, Collinsella_stercoris, Enorma_massiliensis, Adlercreutzia_equolifaciens, Asaccharobacter_celatus, Denitrobacterium_detoxificans, Enterorhabdus_caecimuris, Slackia_exigua, Rubrobacter_radiotolerans, Patulibacter_medicamentivorans, Porphyromonas_endodontalis, Porphyromonas_somerae, Prevotella_histicola, Prevotella_melaninogenica, Prevotella_oris, Capnocytophaga_leadbetteri, Kouleothrix_aurantiaca, Hydrogenibacillus_schlegelii, Gemella_bergeri, Gemella_morbillorum, Brochothrix_campestris, Staphylococcus_schweitzeri, Abiotrophia_defectiva, Abiotrophia_sp_HMSC24B09, Granulicatella_adiacens, Granulicatella_elegans, Enterococcus_faecalis, Lactobacillus_rhamnosus, Streptococcus_anginosus_group, Streptococcus_australis, Streptococcus_cristatus, Streptococcus_gordonii, Streptococcus_infantis, Streptococcus_milleri, Streptococcus_mitis, Streptococcus_parasanguinis, Streptococcus_peroris, Streptococcus_pneumoniae, Streptococcus_pseudopneumoniae, Streptococcus_salivarius, Streptococcus_sanguinis, Streptococcus_sinensis, Streptococcus_sp_A12, Streptococcus_sp_F0442, Streptococcus_sp_HMSC034E03, Streptococcus_sp_HMSC067H01, Streptococcus_sp_HMSC071D03, Streptococcus_sp_HPH0090, Streptococcus_sp_M334, Streptococcus_vestibularis, Eubacterium_brachy, Eubacterium_infirmum, Eubacterium_nodatum, Eubacterium_sulci, Mogibacterium_diversum, Mogibacterium_pumilum, Mogibacterium_timidum, Lachnoanaerobaculum_saburreum, Lachnospiraceae_bacterium_oral_taxon_096, Oribacterium_parvum, Oribacterium_sinus, Oribacterium_sp_oral_taxon_078, Shuttleworthia_satelles, Stomatobaculum_longum, Peptostreptococcus_stomatis, Bulleidia_extructa, Solobacterium_moorei, Selenomonas_noxia, Selenomonas_sputigena, Veillonella_atypica, Veillonella_dispar, Veillonella_infantium, Veillonella_parvula, Veillonella_sp_T11011_6, Anaerococcus_octavius, Parvimonas_micra, Parvimonas_sp_oral_taxon_110, Parvimonas_sp_oral_taxon_393, Peptoniphilus_lacrimalis, Fusobacterium_nucleatum, Leptotrichia_wadei, Gemmata_obscuriglobus, Paludisphaera_borealis, Achromobacter_denitrificans, Achromobacter_ruhlandii, Achromobacter_xylosoxidans, Lautropia_mirabilis, Leptothrix_ochracea, Anaerobiospirillum_thomasii, Cardiobacterium_hominis, Cardiobacterium_valvarum, Alkalilimnicola_ehrlichii, Escherichia_coli, Klebsiella_michiganensis, Klebsiella_oxytoca, Serratia_liquefaciens, Thiohalorhabdus_denitrificans, Pseudomonas_aeruginosa_group, Fretibacterium_fastidiosum, Acholeplasma_oculi, Aspergillus_eucalypticola, Aspergillus_kawachii, Aspergillus_lacticoffeatus, Aspergillus_niger, Aspergillus_phoenicis, Aspergillus_sydowii, Aspergillus_thermomutatus, Aspergillus_tubingensis, Aspergillus_turcosus, Aspergillus_vadensis, Aspergillus_welwitschiae, Candida_albicans, Candida_dubliniensis, Candida_parapsilosis, Candida_tropicalis, Malassezia_globosa, Actinomyces_cardiffensis, Actinomyces_denticolens, Actinomyces_radingae, Actinomyces_turicensis, Varibaculum_cambriense, Bifidobacterium_breve, Gardnerella_vaginalis, Scardovia_inopinata, Corynebacterium_pyruviciproducens, Gordonia_polyisoprenivorans, Mycolicibacterium_fortuitum, Mycolicibacterium_rufum, Nakamurella_silvestris, Atopobium_minutum, Gordonibacter_pamelaeae, Bacteroidetes_oral_taxon_274, Porphyromonas_asaccharolytica, Porphyromonas_canoris, Porphyromonas_catoniae, Porphyromonas_uenonis, Alloprevotella_rava, Alloprevotella_tannerae, Prevotella_buccae, Prevotella_denticola, Prevotella_intermedia, Prevotella_jejuni, Prevotella_nigrescens, Prevotella_oulorum, Prevotella_pallens, Prevotella_pleuritidis, Prevotella_salivae, Prevotella_scopos, Prevotella_sp_F0091, Prevotella_sp_oral_taxon_306, Prevotella_veroralis, Tannerella_forsythia, Tannerella_sp_oral_taxon_808, Tannerella_sp_oral_taxon_HOT_286, Capnocytophaga_gingivalis, Capnocytophaga_granulosa, Capnocytophaga_sputigena, Litorilinea_aerophila, Bacillus_cereus_group, Bacillus_ginsengihumi, Listeria_floridensis, Agitococcus_lubricus, Lactobacillus_fermentum, Lactobacillus_gasseri, Lactobacillus_oris, Lactobacillus_paragasseri, Lactobacillus_reuteri, Lactobacillus_salivarius, Sharpea_azabuensis, Streptococcus_massiliensis, Streptococcus_mutans, Streptococcus_sp_SK643, Streptococcus_sp_oral_taxon_056, Mageeibacillus_indolicus, Oribacterium_asaccharolyticum, Peptostreptococcus_sp_MV1, Megasphaera_micronuciformis, Leptotrichia_sp_oral_taxon_212, Neisseria_flavescens, Neisseria_subflava, Desulfobulbus_oralis, Moraxella_catarrhalis, Stenotrophomonas_maltophilia, Stenotrophomonas_pavanii, Stenotrophomonas_rhizophila, Treponema_socranskii, Bacillus_intestinalis, Listeria_innocua, Listeria_monocytogenes, Staphylococcus_haemolyticus, Enterococcus_avium, Streptococcus_downei, Streptococcus_pyogenes, Streptococcus_salivarius_CAG_79, Streptococcus_sobrinus, Streptococcus_sp_DD11, Streptococcus_sp_HMSC070B10, Streptococcus_sp_NLAE_zl_C503, Streptococcus_thermophilus, Streptococcus_viridans, Eubacterium_saphenum, Eubacterium_rectale, Filifactor_alocis, Eggerthia_catenaformis, Selenomonas_flueggei, Selenomonas_sp_oral_taxon_892, Selenomonas_sp_oral_taxon_920, Dialister_invisus, Dialister_micraerophilus, Negativicoccus_succinicivorans, Veillonella_tobetsuensis, Fusobacterium_periodonticum, Fusobacterium_sp_oral_taxon_370, Leptotrichia_buccalis, Leptotrichia_hofstadii, Leptotrichia_sp_oral_taxon_215, Leptotrichia_sp_oral_taxon_225, Leptotrichia_sp_oral_taxon_498, Leptotrichia_sp_oral_taxon_879, Fimbriiglobus_ruber, Rickettsia_felis, Burkholderia_multivorans, Eikenella_corrodens, Neisseria_elongata, Campylobacter_showae, Haemophilus_parainfluenzae, Candida_orthopsilosis, Rickettsia_typhi, Neisseria_macacae, Neisseria_mucosa, Neisseria_sicca, Nannocystis_exedens, Campylobacter_concisus, Campylobacter_gracilis, Campylobacter_mucosalis, Franconibacter_helveticus, Klebsiella_pneumoniae, Salmonella_enterica, Superficieibacter_electus, Serratia_marcescens, Haemophilus_sp_HMSC71H05, Pseudomonas_fluorescens_group, Pseudomonas_geniculata, Aspergillus_fumigatus, Saccharomyces_cerevisiae, Saccharomyces_cerevisiae_x_Saccharomyces_kudriavzevii, Saccharomyces_kudriavzevii, Cryptococcus_gattii_VGI, Cryptococcus_gattii_VGII, Cryptococcus_gattii_VGIII, Cryptococcus_neoformans
## 2023-03-22 00:45:29 INFO::Total filtered features with variance filtering: 0
## 2023-03-22 00:45:29 INFO::Filtered feature names from variance filtering:
## 2023-03-22 00:45:29 INFO::Running selected normalization method: TSS
## 2023-03-22 00:45:29 INFO::Applying z-score to standardize continuous metadata
## 2023-03-22 00:45:29 INFO::Running selected transform method: LOG
## 2023-03-22 00:45:29 INFO::Running selected analysis method: LM
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 1, Aeriscardovia_aeriphila
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 2, Corynebacterium_accolens
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 3, Corynebacterium_atypicum
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 4, Corynebacterium_aurimucosum
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 5, Corynebacterium_pseudodiphtheriticum
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 6, Corynebacterium_pseudogenitalium
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 7, Rothia_mucilaginosa
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 8, Cutibacterium_acnes
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 9, Cutibacterium_avidum
## 2023-03-22 00:45:29 INFO::Fitting model to feature number 10, Cutibacterium_granulosum
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 11, Propionibacterium_humerusii
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 12, Propionibacterium_namnetense
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 13, Collinsella_aerofaciens
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 14, Collinsella_intestinalis
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 15, X.Collinsella._massiliensis
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 16, Slackia_isoflavoniconvertens
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 17, Thermoleophilum_album
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 18, Gemella_asaccharolytica
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 19, Gemella_haemolysans
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 20, Gemella_sanguinis
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 21, Staphylococcus_argenteus
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 22, Staphylococcus_aureus
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 23, Staphylococcus_capitis
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 24, Staphylococcus_epidermidis
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 25, Dolosigranulum_pigrum
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 26, Streptococcus_oralis
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 27, Limnochorda_pilosa
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 28, Anaerococcus_nagyae
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 29, Finegoldia_magna
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 30, Cupriavidus_sp
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 31, Sutterella_parvirubra
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 32, Malassezia_restricta
## 2023-03-22 00:45:30 INFO::Counting total values for each feature
## 2023-03-22 00:45:30 WARNING::Deleting existing residuals file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/residuals.rds
## 2023-03-22 00:45:30 INFO::Writing residuals to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/residuals.rds
## 2023-03-22 00:45:30 WARNING::Deleting existing fitted file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/fitted.rds
## 2023-03-22 00:45:30 INFO::Writing fitted values to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/fitted.rds
## 2023-03-22 00:45:30 WARNING::Deleting existing ranef file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/ranef.rds
## 2023-03-22 00:45:30 INFO::Writing extracted random effects to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/ranef.rds
## 2023-03-22 00:45:30 INFO::Writing all results to file (ordered by increasing q-values): /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/all_results.tsv
## 2023-03-22 00:45:30 INFO::Writing the significant results (those which are less than or equal to the threshold of 0.250000 ) to file (ordered by increasing q-values): /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/significant_results.tsv
#Checking number of bugs for the main text of the manuscript
cat("Number of differentially abundant bugs by each metadata in NS")
## Number of differentially abundant bugs by each metadata in NS
fit_data_ns$results %>% subset(., .$qval < 0.1) %>% .$metadata %>% table()
## .
## benzonase host_zero log10.Final_reads lypma
## 2 6 9 4
## molysis qiaamp
## 4 1
#bal
# # y ~ log(final reads) + sample_type + treatment
fit_data_bal = Maaslin2(input_data = otu_table(subset_samples(phyloseq_rel_nz, sample_type == "BAL")) %>% t %>% data.frame(),
input_metadata = sample_data(subset_samples(phyloseq_rel_nz, sample_type == "BAL")) %>% data.frame(),
output = "/Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output",
fixed_effects = c("log10.Final_reads","lypma", "benzonase", "host_zero", "molysis", "qiaamp"),
transform = "LOG", #default
normalization = "TSS", # default
random_effects = c("subject_id"),
plot_heatmap = F,
plot_scatter = F)
## [1] "Warning: Deleting existing log file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/maaslin2.log"
## 2023-03-22 00:45:30 INFO::Writing function arguments to log file
## 2023-03-22 00:45:30 INFO::Verifying options selected are valid
## 2023-03-22 00:45:30 INFO::Determining format of input files
## 2023-03-22 00:45:30 INFO::Input format is data samples as rows and metadata samples as rows
## 2023-03-22 00:45:30 INFO::Formula for random effects: expr ~ (1 | subject_id)
## 2023-03-22 00:45:30 INFO::Formula for fixed effects: expr ~ log10.Final_reads + lypma + benzonase + host_zero + molysis + qiaamp
## 2023-03-22 00:45:30 INFO::Filter data based on min abundance and min prevalence
## 2023-03-22 00:45:30 INFO::Total samples in data: 28
## 2023-03-22 00:45:30 INFO::Min samples required with min abundance for a feature not to be filtered: 2.800000
## 2023-03-22 00:45:30 INFO::Total filtered features: 261
## 2023-03-22 00:45:30 INFO::Filtered feature names from abundance and prevalence filtering: Actinobaculum_sp_oral_taxon_183, Actinomyces_georgiae, Actinomyces_graevenitzii, Actinomyces_hongkongensis, Actinomyces_johnsonii, Actinomyces_massiliensis, Actinomyces_meyeri, Actinomyces_oris, Actinomyces_sp_HMSC035G02, Actinomyces_sp_ICM47, Actinomyces_sp_oral_taxon_170, Actinomyces_sp_oral_taxon_180, Actinomyces_sp_oral_taxon_414, Actinomyces_sp_oral_taxon_448, Actinomyces_sp_oral_taxon_897, Actinomyces_viscosus, Aeriscardovia_aeriphila, Alloscardovia_omnicolens, Bifidobacterium_dentium, Bifidobacterium_longum, Parascardovia_denticolens, Corynebacterium_accolens, Corynebacterium_afermentans, Corynebacterium_atypicum, Corynebacterium_aurimucosum, Corynebacterium_durum, Corynebacterium_kroppenstedtii, Corynebacterium_matruchotii, Corynebacterium_pseudodiphtheriticum, Corynebacterium_pseudogenitalium, Mycobacterium_chimaera, Mycobacterium_intracellulare, Micrococcus_luteus, Rothia_aeria, Rothia_dentocariosa, Tropheryma_whipplei, Cutibacterium_avidum, Cutibacterium_granulosum, Propionibacterium_humerusii, Propionibacterium_namnetense, Propionibacterium_acidifaciens, Pseudopropionibacterium_propionicum, Olsenella_profusa, Olsenella_uli, Collinsella_aerofaciens, Collinsella_intestinalis, Collinsella_stercoris, Enorma_massiliensis, X.Collinsella._massiliensis, Adlercreutzia_equolifaciens, Asaccharobacter_celatus, Denitrobacterium_detoxificans, Enterorhabdus_caecimuris, Slackia_exigua, Slackia_isoflavoniconvertens, Rubrobacter_radiotolerans, Patulibacter_medicamentivorans, Thermoleophilum_album, Porphyromonas_endodontalis, Porphyromonas_somerae, Prevotella_histicola, Capnocytophaga_leadbetteri, Kouleothrix_aurantiaca, Hydrogenibacillus_schlegelii, Brochothrix_campestris, Staphylococcus_capitis, Abiotrophia_defectiva, Abiotrophia_sp_HMSC24B09, Dolosigranulum_pigrum, Granulicatella_adiacens, Granulicatella_elegans, Lactobacillus_rhamnosus, Streptococcus_australis, Streptococcus_cristatus, Streptococcus_gordonii, Streptococcus_infantis, Streptococcus_milleri, Streptococcus_peroris, Streptococcus_pseudopneumoniae, Streptococcus_salivarius, Streptococcus_sanguinis, Streptococcus_sinensis, Streptococcus_sp_A12, Streptococcus_sp_F0442, Streptococcus_sp_HMSC034E03, Streptococcus_sp_HMSC067H01, Streptococcus_sp_HMSC071D03, Streptococcus_sp_HPH0090, Streptococcus_sp_M334, Streptococcus_vestibularis, Eubacterium_brachy, Eubacterium_nodatum, Eubacterium_sulci, Mogibacterium_diversum, Mogibacterium_pumilum, Mogibacterium_timidum, Lachnoanaerobaculum_saburreum, Lachnospiraceae_bacterium_oral_taxon_096, Oribacterium_parvum, Oribacterium_sinus, Oribacterium_sp_oral_taxon_078, Shuttleworthia_satelles, Stomatobaculum_longum, Solobacterium_moorei, Limnochorda_pilosa, Selenomonas_noxia, Selenomonas_sputigena, Veillonella_atypica, Veillonella_dispar, Veillonella_infantium, Veillonella_sp_T11011_6, Anaerococcus_nagyae, Anaerococcus_octavius, Parvimonas_sp_oral_taxon_110, Parvimonas_sp_oral_taxon_393, Peptoniphilus_lacrimalis, Leptotrichia_wadei, Gemmata_obscuriglobus, Paludisphaera_borealis, Achromobacter_denitrificans, Achromobacter_ruhlandii, Achromobacter_xylosoxidans, Lautropia_mirabilis, Leptothrix_ochracea, Anaerobiospirillum_thomasii, Cardiobacterium_hominis, Alkalilimnicola_ehrlichii, Escherichia_coli, Klebsiella_michiganensis, Klebsiella_oxytoca, Serratia_liquefaciens, Thiohalorhabdus_denitrificans, Pseudomonas_aeruginosa_group, Fretibacterium_fastidiosum, Acholeplasma_oculi, Aspergillus_eucalypticola, Aspergillus_kawachii, Aspergillus_lacticoffeatus, Aspergillus_niger, Aspergillus_phoenicis, Aspergillus_sydowii, Aspergillus_thermomutatus, Aspergillus_tubingensis, Aspergillus_turcosus, Aspergillus_vadensis, Aspergillus_welwitschiae, Candida_tropicalis, Malassezia_globosa, Malassezia_restricta, Actinomyces_cardiffensis, Actinomyces_denticolens, Actinomyces_turicensis, Bifidobacterium_breve, Gardnerella_vaginalis, Scardovia_inopinata, Corynebacterium_pyruviciproducens, Gordonia_polyisoprenivorans, Mycolicibacterium_fortuitum, Mycolicibacterium_rufum, Nakamurella_silvestris, Atopobium_minutum, Gordonibacter_pamelaeae, Bacteroidetes_oral_taxon_274, Porphyromonas_canoris, Porphyromonas_catoniae, Porphyromonas_uenonis, Prevotella_denticola, Prevotella_intermedia, Prevotella_nigrescens, Prevotella_oulorum, Prevotella_pallens, Prevotella_pleuritidis, Prevotella_scopos, Prevotella_sp_F0091, Tannerella_forsythia, Tannerella_sp_oral_taxon_808, Tannerella_sp_oral_taxon_HOT_286, Capnocytophaga_gingivalis, Capnocytophaga_granulosa, Capnocytophaga_sputigena, Litorilinea_aerophila, Bacillus_cereus_group, Bacillus_ginsengihumi, Listeria_floridensis, Agitococcus_lubricus, Lactobacillus_gasseri, Lactobacillus_oris, Lactobacillus_paragasseri, Lactobacillus_reuteri, Lactobacillus_salivarius, Sharpea_azabuensis, Streptococcus_massiliensis, Streptococcus_mutans, Streptococcus_sp_SK643, Streptococcus_sp_oral_taxon_056, Mageeibacillus_indolicus, Oribacterium_asaccharolyticum, Leptotrichia_sp_oral_taxon_212, Neisseria_subflava, Desulfobulbus_oralis, Moraxella_catarrhalis, Stenotrophomonas_maltophilia, Stenotrophomonas_pavanii, Stenotrophomonas_rhizophila, Treponema_socranskii, Bacillus_intestinalis, Listeria_innocua, Listeria_monocytogenes, Streptococcus_downei, Streptococcus_salivarius_CAG_79, Streptococcus_sobrinus, Streptococcus_sp_DD11, Streptococcus_sp_HMSC070B10, Streptococcus_sp_NLAE_zl_C503, Streptococcus_viridans, Eubacterium_saphenum, Eubacterium_rectale, Filifactor_alocis, Eggerthia_catenaformis, Selenomonas_flueggei, Selenomonas_sp_oral_taxon_892, Selenomonas_sp_oral_taxon_920, Dialister_micraerophilus, Veillonella_tobetsuensis, Fusobacterium_periodonticum, Fusobacterium_sp_oral_taxon_370, Leptotrichia_buccalis, Leptotrichia_hofstadii, Leptotrichia_sp_oral_taxon_215, Leptotrichia_sp_oral_taxon_225, Leptotrichia_sp_oral_taxon_498, Leptotrichia_sp_oral_taxon_879, Fimbriiglobus_ruber, Rickettsia_felis, Burkholderia_multivorans, Eikenella_corrodens, Neisseria_elongata, Campylobacter_showae, Haemophilus_parainfluenzae, Rickettsia_typhi, Neisseria_macacae, Neisseria_mucosa, Neisseria_sicca, Nannocystis_exedens, Campylobacter_gracilis, Campylobacter_mucosalis, Franconibacter_helveticus, Klebsiella_pneumoniae, Salmonella_enterica, Superficieibacter_electus, Haemophilus_sp_HMSC71H05, Pseudomonas_fluorescens_group, Pseudomonas_geniculata, Aspergillus_fumigatus, Saccharomyces_cerevisiae, Saccharomyces_cerevisiae_x_Saccharomyces_kudriavzevii, Saccharomyces_kudriavzevii, Cryptococcus_gattii_VGI, Cryptococcus_gattii_VGII, Cryptococcus_gattii_VGIII, Cryptococcus_neoformans
## 2023-03-22 00:45:30 INFO::Total filtered features with variance filtering: 0
## 2023-03-22 00:45:30 INFO::Filtered feature names from variance filtering:
## 2023-03-22 00:45:30 INFO::Running selected normalization method: TSS
## 2023-03-22 00:45:30 INFO::Applying z-score to standardize continuous metadata
## 2023-03-22 00:45:30 INFO::Running selected transform method: LOG
## 2023-03-22 00:45:30 INFO::Running selected analysis method: LM
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 1, Actinomyces_naeslundii
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 2, Actinomyces_odontolyticus
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 3, Actinomyces_sp_HPA0247
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 4, Actinomyces_sp_S6_Spd3
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 5, Actinomyces_sp_oral_taxon_181
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 6, Scardovia_wiggsiae
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 7, Corynebacterium_tuberculostearicum
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 8, Rothia_mucilaginosa
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 9, Cutibacterium_acnes
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 10, Atopobium_deltae
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 11, Atopobium_parvulum
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 12, Atopobium_rimae
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 13, Olsenella_scatoligenes
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 14, Prevotella_melaninogenica
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 15, Prevotella_oris
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 16, Gemella_asaccharolytica
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 17, Gemella_bergeri
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 18, Gemella_haemolysans
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 19, Gemella_morbillorum
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 20, Gemella_sanguinis
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 21, Staphylococcus_argenteus
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 22, Staphylococcus_aureus
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 23, Staphylococcus_epidermidis
## 2023-03-22 00:45:30 INFO::Fitting model to feature number 24, Staphylococcus_schweitzeri
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 25, Enterococcus_faecalis
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 26, Streptococcus_anginosus_group
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 27, Streptococcus_mitis
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 28, Streptococcus_oralis
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 29, Streptococcus_parasanguinis
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 30, Streptococcus_pneumoniae
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 31, Eubacterium_infirmum
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 32, Peptostreptococcus_stomatis
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 33, Bulleidia_extructa
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 34, Veillonella_parvula
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 35, Finegoldia_magna
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 36, Parvimonas_micra
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 37, Fusobacterium_nucleatum
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 38, Cupriavidus_sp
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 39, Sutterella_parvirubra
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 40, Cardiobacterium_valvarum
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 41, Candida_albicans
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 42, Candida_dubliniensis
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 43, Candida_parapsilosis
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 44, Actinomyces_radingae
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 45, Varibaculum_cambriense
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 46, Porphyromonas_asaccharolytica
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 47, Alloprevotella_rava
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 48, Alloprevotella_tannerae
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 49, Prevotella_buccae
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 50, Prevotella_jejuni
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 51, Prevotella_salivae
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 52, Prevotella_sp_oral_taxon_306
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 53, Prevotella_veroralis
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 54, Lactobacillus_fermentum
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 55, Peptostreptococcus_sp_MV1
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 56, Megasphaera_micronuciformis
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 57, Neisseria_flavescens
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 58, Staphylococcus_haemolyticus
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 59, Enterococcus_avium
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 60, Streptococcus_pyogenes
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 61, Streptococcus_thermophilus
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 62, Dialister_invisus
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 63, Negativicoccus_succinicivorans
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 64, Candida_orthopsilosis
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 65, Campylobacter_concisus
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 66, Serratia_marcescens
## 2023-03-22 00:45:31 INFO::Counting total values for each feature
## 2023-03-22 00:45:31 WARNING::Deleting existing residuals file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/residuals.rds
## 2023-03-22 00:45:31 INFO::Writing residuals to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/residuals.rds
## 2023-03-22 00:45:31 WARNING::Deleting existing fitted file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/fitted.rds
## 2023-03-22 00:45:31 INFO::Writing fitted values to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/fitted.rds
## 2023-03-22 00:45:31 WARNING::Deleting existing ranef file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/ranef.rds
## 2023-03-22 00:45:31 INFO::Writing extracted random effects to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/ranef.rds
## 2023-03-22 00:45:31 INFO::Writing all results to file (ordered by increasing q-values): /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/all_results.tsv
## 2023-03-22 00:45:31 INFO::Writing the significant results (those which are less than or equal to the threshold of 0.250000 ) to file (ordered by increasing q-values): /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/significant_results.tsv
cat("Number of differentially abundant bugs by each metadata in BAL")
## Number of differentially abundant bugs by each metadata in BAL
fit_data_bal$results %>% subset(., .$qval < 0.1) %>% .$metadata %>% table()
## < table of extent 0 >
#sputum
# # y ~ log(final reads) + sample_type + treatment
fit_data_spt = Maaslin2(input_data = otu_table(subset_samples(phyloseq_rel_nz, sample_type == "Sputum")) %>% t %>% data.frame(),
input_metadata = sample_data(subset_samples(phyloseq_rel_nz, sample_type == "Sputum")) %>% data.frame(),
output = "/Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output",
fixed_effects = c("log10.Final_reads","lypma", "benzonase", "host_zero", "molysis", "qiaamp"),
transform = "LOG", #default
normalization = "TSS", # default
random_effects = c("subject_id"),
plot_heatmap = F,
plot_scatter = F)
## [1] "Warning: Deleting existing log file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/maaslin2.log"
## 2023-03-22 00:45:31 INFO::Writing function arguments to log file
## 2023-03-22 00:45:31 INFO::Verifying options selected are valid
## 2023-03-22 00:45:31 INFO::Determining format of input files
## 2023-03-22 00:45:31 INFO::Input format is data samples as rows and metadata samples as rows
## 2023-03-22 00:45:31 INFO::Formula for random effects: expr ~ (1 | subject_id)
## 2023-03-22 00:45:31 INFO::Formula for fixed effects: expr ~ log10.Final_reads + lypma + benzonase + host_zero + molysis + qiaamp
## 2023-03-22 00:45:31 INFO::Filter data based on min abundance and min prevalence
## 2023-03-22 00:45:31 INFO::Total samples in data: 30
## 2023-03-22 00:45:31 INFO::Min samples required with min abundance for a feature not to be filtered: 3.000000
## 2023-03-22 00:45:31 INFO::Total filtered features: 155
## 2023-03-22 00:45:31 INFO::Filtered feature names from abundance and prevalence filtering: Parascardovia_denticolens, Corynebacterium_accolens, Corynebacterium_afermentans, Corynebacterium_aurimucosum, Corynebacterium_kroppenstedtii, Corynebacterium_pseudodiphtheriticum, Corynebacterium_pseudogenitalium, Corynebacterium_tuberculostearicum, Mycobacterium_chimaera, Mycobacterium_intracellulare, Micrococcus_luteus, Tropheryma_whipplei, Cutibacterium_avidum, Cutibacterium_granulosum, Propionibacterium_humerusii, Propionibacterium_namnetense, Patulibacter_medicamentivorans, Prevotella_oris, Gemella_asaccharolytica, Staphylococcus_capitis, Dolosigranulum_pigrum, Streptococcus_sinensis, Selenomonas_noxia, Selenomonas_sputigena, Anaerococcus_nagyae, Anaerococcus_octavius, Finegoldia_magna, Lautropia_mirabilis, Leptothrix_ochracea, Klebsiella_michiganensis, Klebsiella_oxytoca, Serratia_liquefaciens, Aspergillus_eucalypticola, Aspergillus_kawachii, Aspergillus_lacticoffeatus, Aspergillus_niger, Aspergillus_phoenicis, Aspergillus_sydowii, Aspergillus_thermomutatus, Aspergillus_tubingensis, Aspergillus_turcosus, Aspergillus_vadensis, Aspergillus_welwitschiae, Candida_albicans, Candida_dubliniensis, Candida_tropicalis, Malassezia_globosa, Malassezia_restricta, Actinomyces_cardiffensis, Actinomyces_denticolens, Actinomyces_radingae, Actinomyces_turicensis, Varibaculum_cambriense, Gardnerella_vaginalis, Scardovia_inopinata, Corynebacterium_pyruviciproducens, Gordonia_polyisoprenivorans, Mycolicibacterium_fortuitum, Mycolicibacterium_rufum, Atopobium_minutum, Gordonibacter_pamelaeae, Bacteroidetes_oral_taxon_274, Porphyromonas_asaccharolytica, Porphyromonas_canoris, Porphyromonas_uenonis, Alloprevotella_rava, Alloprevotella_tannerae, Prevotella_buccae, Prevotella_denticola, Prevotella_intermedia, Prevotella_nigrescens, Prevotella_oulorum, Prevotella_pleuritidis, Prevotella_scopos, Prevotella_sp_F0091, Prevotella_sp_oral_taxon_306, Prevotella_veroralis, Tannerella_sp_oral_taxon_808, Capnocytophaga_granulosa, Capnocytophaga_sputigena, Bacillus_cereus_group, Bacillus_ginsengihumi, Listeria_floridensis, Agitococcus_lubricus, Lactobacillus_fermentum, Lactobacillus_oris, Lactobacillus_paragasseri, Lactobacillus_reuteri, Lactobacillus_salivarius, Streptococcus_massiliensis, Streptococcus_sp_SK643, Mageeibacillus_indolicus, Peptostreptococcus_sp_MV1, Leptotrichia_sp_oral_taxon_212, Desulfobulbus_oralis, Moraxella_catarrhalis, Treponema_socranskii, Bacillus_intestinalis, Listeria_innocua, Listeria_monocytogenes, Staphylococcus_haemolyticus, Enterococcus_avium, Streptococcus_downei, Streptococcus_pyogenes, Streptococcus_salivarius_CAG_79, Streptococcus_sp_DD11, Streptococcus_sp_HMSC070B10, Streptococcus_sp_NLAE_zl_C503, Streptococcus_thermophilus, Eubacterium_rectale, Selenomonas_flueggei, Selenomonas_sp_oral_taxon_892, Selenomonas_sp_oral_taxon_920, Dialister_invisus, Dialister_micraerophilus, Negativicoccus_succinicivorans, Veillonella_tobetsuensis, Fusobacterium_periodonticum, Fusobacterium_sp_oral_taxon_370, Leptotrichia_buccalis, Leptotrichia_hofstadii, Leptotrichia_sp_oral_taxon_225, Leptotrichia_sp_oral_taxon_498, Leptotrichia_sp_oral_taxon_879, Fimbriiglobus_ruber, Rickettsia_felis, Burkholderia_multivorans, Eikenella_corrodens, Neisseria_elongata, Campylobacter_showae, Candida_orthopsilosis, Rickettsia_typhi, Neisseria_macacae, Neisseria_mucosa, Neisseria_sicca, Nannocystis_exedens, Campylobacter_concisus, Campylobacter_gracilis, Campylobacter_mucosalis, Franconibacter_helveticus, Klebsiella_pneumoniae, Salmonella_enterica, Superficieibacter_electus, Serratia_marcescens, Haemophilus_sp_HMSC71H05, Pseudomonas_fluorescens_group, Pseudomonas_geniculata, Aspergillus_fumigatus, Saccharomyces_cerevisiae, Saccharomyces_cerevisiae_x_Saccharomyces_kudriavzevii, Saccharomyces_kudriavzevii, Cryptococcus_gattii_VGI, Cryptococcus_gattii_VGII, Cryptococcus_gattii_VGIII, Cryptococcus_neoformans
## 2023-03-22 00:45:31 INFO::Total filtered features with variance filtering: 0
## 2023-03-22 00:45:31 INFO::Filtered feature names from variance filtering:
## 2023-03-22 00:45:31 INFO::Running selected normalization method: TSS
## 2023-03-22 00:45:31 INFO::Applying z-score to standardize continuous metadata
## 2023-03-22 00:45:31 INFO::Running selected transform method: LOG
## 2023-03-22 00:45:31 INFO::Running selected analysis method: LM
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 1, Actinobaculum_sp_oral_taxon_183
## 2023-03-22 00:45:31 INFO::Fitting model to feature number 2, Actinomyces_georgiae
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 3, Actinomyces_graevenitzii
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 4, Actinomyces_hongkongensis
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 5, Actinomyces_johnsonii
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 6, Actinomyces_massiliensis
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 7, Actinomyces_meyeri
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 8, Actinomyces_naeslundii
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 9, Actinomyces_odontolyticus
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 10, Actinomyces_oris
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 11, Actinomyces_sp_HMSC035G02
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 12, Actinomyces_sp_HPA0247
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 13, Actinomyces_sp_ICM47
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 14, Actinomyces_sp_S6_Spd3
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 15, Actinomyces_sp_oral_taxon_170
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 16, Actinomyces_sp_oral_taxon_180
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 17, Actinomyces_sp_oral_taxon_181
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 18, Actinomyces_sp_oral_taxon_414
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 19, Actinomyces_sp_oral_taxon_448
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 20, Actinomyces_sp_oral_taxon_897
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 21, Actinomyces_viscosus
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 22, Aeriscardovia_aeriphila
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 23, Alloscardovia_omnicolens
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 24, Bifidobacterium_dentium
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 25, Bifidobacterium_longum
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 26, Scardovia_wiggsiae
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 27, Corynebacterium_atypicum
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 28, Corynebacterium_durum
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 29, Corynebacterium_matruchotii
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 30, Rothia_aeria
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 31, Rothia_dentocariosa
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 32, Rothia_mucilaginosa
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 33, Cutibacterium_acnes
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 34, Propionibacterium_acidifaciens
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 35, Pseudopropionibacterium_propionicum
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 36, Atopobium_deltae
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 37, Atopobium_parvulum
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 38, Atopobium_rimae
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 39, Olsenella_profusa
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 40, Olsenella_scatoligenes
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 41, Olsenella_uli
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 42, Collinsella_aerofaciens
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 43, Collinsella_intestinalis
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 44, Collinsella_stercoris
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 45, Enorma_massiliensis
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 46, X.Collinsella._massiliensis
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 47, Adlercreutzia_equolifaciens
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 48, Asaccharobacter_celatus
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 49, Denitrobacterium_detoxificans
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 50, Enterorhabdus_caecimuris
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 51, Slackia_exigua
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 52, Slackia_isoflavoniconvertens
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 53, Rubrobacter_radiotolerans
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 54, Thermoleophilum_album
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 55, Porphyromonas_endodontalis
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 56, Porphyromonas_somerae
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 57, Prevotella_histicola
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 58, Prevotella_melaninogenica
## 2023-03-22 00:45:32 INFO::Fitting model to feature number 59, Capnocytophaga_leadbetteri
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 60, Kouleothrix_aurantiaca
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 61, Hydrogenibacillus_schlegelii
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 62, Gemella_bergeri
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 63, Gemella_haemolysans
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 64, Gemella_morbillorum
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 65, Gemella_sanguinis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 66, Brochothrix_campestris
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 67, Staphylococcus_argenteus
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 68, Staphylococcus_aureus
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 69, Staphylococcus_epidermidis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 70, Staphylococcus_schweitzeri
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 71, Abiotrophia_defectiva
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 72, Abiotrophia_sp_HMSC24B09
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 73, Granulicatella_adiacens
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 74, Granulicatella_elegans
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 75, Enterococcus_faecalis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 76, Lactobacillus_rhamnosus
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 77, Streptococcus_anginosus_group
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 78, Streptococcus_australis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 79, Streptococcus_cristatus
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 80, Streptococcus_gordonii
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 81, Streptococcus_infantis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 82, Streptococcus_milleri
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 83, Streptococcus_mitis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 84, Streptococcus_oralis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 85, Streptococcus_parasanguinis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 86, Streptococcus_peroris
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 87, Streptococcus_pneumoniae
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 88, Streptococcus_pseudopneumoniae
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 89, Streptococcus_salivarius
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 90, Streptococcus_sanguinis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 91, Streptococcus_sp_A12
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 92, Streptococcus_sp_F0442
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 93, Streptococcus_sp_HMSC034E03
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 94, Streptococcus_sp_HMSC067H01
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 95, Streptococcus_sp_HMSC071D03
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 96, Streptococcus_sp_HPH0090
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 97, Streptococcus_sp_M334
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 98, Streptococcus_vestibularis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 99, Eubacterium_brachy
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 100, Eubacterium_infirmum
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 101, Eubacterium_nodatum
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 102, Eubacterium_sulci
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 103, Mogibacterium_diversum
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 104, Mogibacterium_pumilum
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 105, Mogibacterium_timidum
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 106, Lachnoanaerobaculum_saburreum
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 107, Lachnospiraceae_bacterium_oral_taxon_096
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 108, Oribacterium_parvum
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 109, Oribacterium_sinus
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 110, Oribacterium_sp_oral_taxon_078
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 111, Shuttleworthia_satelles
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 112, Stomatobaculum_longum
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 113, Peptostreptococcus_stomatis
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 114, Bulleidia_extructa
## 2023-03-22 00:45:33 INFO::Fitting model to feature number 115, Solobacterium_moorei
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 116, Limnochorda_pilosa
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 117, Veillonella_atypica
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 118, Veillonella_dispar
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 119, Veillonella_infantium
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 120, Veillonella_parvula
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 121, Veillonella_sp_T11011_6
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 122, Parvimonas_micra
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 123, Parvimonas_sp_oral_taxon_110
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 124, Parvimonas_sp_oral_taxon_393
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 125, Peptoniphilus_lacrimalis
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 126, Fusobacterium_nucleatum
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 127, Leptotrichia_wadei
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 128, Gemmata_obscuriglobus
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 129, Paludisphaera_borealis
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 130, Achromobacter_denitrificans
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 131, Achromobacter_ruhlandii
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 132, Achromobacter_xylosoxidans
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 133, Cupriavidus_sp
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 134, Sutterella_parvirubra
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 135, Anaerobiospirillum_thomasii
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 136, Cardiobacterium_hominis
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 137, Cardiobacterium_valvarum
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 138, Alkalilimnicola_ehrlichii
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 139, Escherichia_coli
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 140, Thiohalorhabdus_denitrificans
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 141, Pseudomonas_aeruginosa_group
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 142, Fretibacterium_fastidiosum
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 143, Acholeplasma_oculi
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 144, Candida_parapsilosis
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 145, Bifidobacterium_breve
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 146, Nakamurella_silvestris
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 147, Porphyromonas_catoniae
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 148, Prevotella_jejuni
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 149, Prevotella_pallens
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 150, Prevotella_salivae
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 151, Tannerella_forsythia
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 152, Tannerella_sp_oral_taxon_HOT_286
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 153, Capnocytophaga_gingivalis
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 154, Litorilinea_aerophila
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 155, Lactobacillus_gasseri
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 156, Sharpea_azabuensis
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 157, Streptococcus_mutans
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 158, Streptococcus_sp_oral_taxon_056
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 159, Oribacterium_asaccharolyticum
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 160, Megasphaera_micronuciformis
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 161, Neisseria_flavescens
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 162, Neisseria_subflava
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 163, Stenotrophomonas_maltophilia
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 164, Stenotrophomonas_pavanii
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 165, Stenotrophomonas_rhizophila
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 166, Streptococcus_sobrinus
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 167, Streptococcus_viridans
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 168, Eubacterium_saphenum
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 169, Filifactor_alocis
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 170, Eggerthia_catenaformis
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 171, Leptotrichia_sp_oral_taxon_215
## 2023-03-22 00:45:34 INFO::Fitting model to feature number 172, Haemophilus_parainfluenzae
## 2023-03-22 00:45:35 INFO::Counting total values for each feature
## 2023-03-22 00:45:35 WARNING::Deleting existing residuals file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/residuals.rds
## 2023-03-22 00:45:35 INFO::Writing residuals to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/residuals.rds
## 2023-03-22 00:45:35 WARNING::Deleting existing fitted file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/fitted.rds
## 2023-03-22 00:45:35 INFO::Writing fitted values to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/fitted.rds
## 2023-03-22 00:45:35 WARNING::Deleting existing ranef file: /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/ranef.rds
## 2023-03-22 00:45:35 INFO::Writing extracted random effects to file /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/ranef.rds
## 2023-03-22 00:45:35 INFO::Writing all results to file (ordered by increasing q-values): /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/all_results.tsv
## 2023-03-22 00:45:35 INFO::Writing the significant results (those which are less than or equal to the threshold of 0.250000 ) to file (ordered by increasing q-values): /Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_host_dna_depletion/Data/maaslin_output/significant_results.tsv
cat("Number of differentially abundant bugs by each metadata in sputum")
## Number of differentially abundant bugs by each metadata in sputum
fit_data_spt$results %>% subset(., .$qval < 0.1) %>% .$metadata %>% table()
## .
## benzonase host_zero log10.Final_reads lypma
## 8 5 1 8
## molysis qiaamp
## 8 7
cat("\nAdding control data to both all and stratified MaAslin will help identify the actual contaminant\n")
##
## Adding control data to both all and stratified MaAslin will help identify the actual contaminant
#Making significance table for figure
# Define a function to make species names italicized
species_italic <- function(data) {
names <- gsub("_", " ", rownames(data))
names <- gsub("[]]|[[]", "", names)
names <- gsub(" sp", " sp.", names)
names <- gsub(" sp.", "* sp.", names)
names <- gsub(" group", "* group.", names)
names <- ifelse(grepl("[*]", names), paste("*", names, sep = ""), paste("*", names, "*", sep = ""))
rownames(data) <- names
data
}
#volcano plot of MaAslin with all data
ggplot(maaslin_all$results, aes(y = -log10(qval), x = coef, col = metadata)) +
theme_classic(base_family = "serif") +
#labs(tag = "A") +
ggtitle("MaAslin with treatment types")+
geom_point(size = 2) +
xlab("MaAslin coefficient") +
ylab("-log<sub>10</sub>(*q*-value)") +
geom_hline(yintercept = 1, col = "gray") +
geom_vline(xintercept = 0, col = "gray") +
geom_richtext(aes( 4, 8, label = "*q*-value = 0.1, fold-change = 0", vjust = -1, fontface = 1), col = "grey", size = 3, family = "serif") +
theme(legend.position = "top", axis.title.y = ggtext::element_markdown()) +
scale_color_manual(values = c("#4daf4a", "#984ea3", "#f781bf", "#377eb8", "#ff7f00", "#ffff33", "#a65628")) +
guides(col = guide_legend(title = "Fixed effects", title.position = "top", nrow = 3))
# Make a significance table for each figure (top 20 taxa)
fit_data_spt$results
make_sig_table <- function(data) {
sig_data <- spread(data$results[order(data$results$qval), c("feature", "metadata", "qval")], metadata, qval)
sig_data$min <- apply(sig_data, 1, FUN = min)
sig_data <- sig_data[order(sig_data$min),] %>% .[1:20,]
sig_data[["feature"]] <- ifelse(sig_data[["feature"]] == "X.Collinsella._massiliensis", "[Collinsella]_massiliensis", sig_data[["feature"]])
sig_data_italic <- sig_data %>% rownames_to_column(var = "-") %>% column_to_rownames(var = "feature") %>% species_italic %>% select(-c("-", "min", "log10.Final_reads"))
sig_data_sig <- ifelse(sig_data_italic < 0.1, "*", NA)
return(list(data = sig_data, data_italic = sig_data_italic, data_sig = sig_data_sig))
}
fit_data_bal <- make_sig_table(fit_data_bal)
fit_data_ns <- make_sig_table(fit_data_ns)
fit_data_spt <- make_sig_table(fit_data_spt)
fit_data_ns$data
cat("\n\nStratified Masslin - significant of each but by treatment in BAL\n")
##
##
## Stratified Masslin - significant of each but by treatment in BAL
fit_data_bal$data_sig
## benzonase host_zero lypma molysis qiaamp
## *Actinomyces radingae* NA NA NA NA NA
## *Bulleidia extructa* NA NA NA NA NA
## *Candida dubliniensis* NA NA NA NA NA
## *Enterococcus avium* NA NA NA NA NA
## *Finegoldia magna* NA NA NA NA NA
## *Gemella morbillorum* NA NA NA NA NA
## *Peptostreptococcus* sp. MV1 NA NA NA NA NA
## *Peptostreptococcus stomatis* NA NA NA NA NA
## *Rothia mucilaginosa* NA NA NA NA NA
## *Staphylococcus epidermidis* NA NA NA NA NA
## *Staphylococcus schweitzeri* NA NA NA NA NA
## *Varibaculum cambriense* NA NA NA NA NA
## *Actinomyces naeslundii* NA NA NA NA NA
## *Corynebacterium tuberculostearicum* NA NA NA NA NA
## *Gemella asaccharolytica* NA NA NA NA NA
## *Gemella haemolysans* NA NA NA NA NA
## *Gemella sanguinis* NA NA NA NA NA
## *Negativicoccus succinicivorans* NA NA NA NA NA
## *Sutterella parvirubra* NA NA NA NA NA
## *Prevotella oris* NA NA NA NA NA
cat("\n\nStratified Masslin - significant of each but by treatment in nasal swab\n")
##
##
## Stratified Masslin - significant of each but by treatment in nasal swab
fit_data_ns$data_sig
## benzonase host_zero lypma molysis qiaamp
## *Cupriavidus* sp. "*" "*" "*" "*" NA
## *Streptococcus oralis* NA NA NA "*" NA
## *Staphylococcus epidermidis* NA NA "*" NA NA
## *Sutterella parvirubra* NA "*" "*" NA NA
## *Aeriscardovia aeriphila* NA NA NA NA NA
## *Rothia mucilaginosa* NA "*" NA "*" NA
## *Corynebacterium atypicum* NA NA NA NA NA
## *Propionibacterium namnetense* NA NA NA NA NA
## *Slackia isoflavoniconvertens* NA NA NA NA NA
## *Cutibacterium avidum* NA NA NA NA NA
## *Staphylococcus aureus* NA NA NA NA NA
## *Staphylococcus argenteus* NA "*" NA NA NA
## *Collinsella massiliensis* NA NA NA NA "*"
## *Limnochorda pilosa* NA NA NA "*" NA
## *Corynebacterium accolens* NA NA "*" NA NA
## *Malassezia restricta* "*" NA NA NA NA
## *Collinsella intestinalis* NA "*" NA NA NA
## *Finegoldia magna* NA "*" NA NA NA
## *Gemella haemolysans* NA NA NA NA NA
## *Thermoleophilum album* NA NA NA NA NA
cat("\n\nStratified Masslin - significant of each but by treatment in Sputum\n")
##
##
## Stratified Masslin - significant of each but by treatment in Sputum
fit_data_spt$data_sig
## benzonase host_zero lypma molysis qiaamp
## *Streptococcus infantis* "*" NA "*" NA NA
## *Gemella haemolysans* "*" NA "*" NA NA
## *Eubacterium sulci* "*" "*" NA "*" "*"
## *Granulicatella elegans* "*" NA "*" NA NA
## *Rothia dentocariosa* "*" NA "*" NA NA
## *Streptococcus* sp. F0442 "*" NA "*" NA NA
## *Bifidobacterium longum* "*" "*" NA "*" "*"
## *Pseudomonas aeruginosa* group. NA NA NA "*" "*"
## *Neisseria flavescens* "*" NA "*" NA NA
## *Rothia mucilaginosa* NA NA "*" NA NA
## *Actinomyces naeslundii* NA "*" NA "*" "*"
## *Denitrobacterium detoxificans* NA "*" NA "*" "*"
## *Actinomyces* sp. oral taxon 897 NA "*" NA "*" "*"
## *Prevotella melaninogenica* NA NA NA "*" "*"
## *Gemmata obscuriglobus* NA NA NA "*" NA
## *Streptococcus gordonii* NA NA "*" NA NA
## *Corynebacterium durum* NA NA NA NA NA
## *Actinomyces* sp. HPA0247 NA NA NA NA NA
## *Streptococcus sanguinis* NA NA NA NA NA
## *Actinomyces massiliensis* NA NA NA NA NA
=======
Results After adding control data, MaAslin needs to be reanalyzed. Adding controls (mock communities) for each treatment group will show more statistically valid results in y ~ log(final reads) + sample_type + treatment, (re = subject_id))
A6. Decontam - stratified by treatment
--> both stratified and nonstratified
input of DNA concentration: 16S qPCR data
https://github.com/benjjneb/decontam/issues/33
Ben Callahan: But in the more limited testing on qPCR data the method still seems to work, and other publications report strong patterns of inverse frequency of contaminants using qPCR data - which is the pattern the frequency method relies on.
Strategy:
2.4.1. run decontam for all samples (common contaminants, by extraction)
2.4.2. stratify decontam analysis per each treatment method (contaminants by depletion methods)
# Decontam package --------------------------------------------------------
# common contaminants across all the treatment methods
#Decontam - were there any contaminants?#
sample_data(phyloseq$phyloseq_rel)$is.neg <- grepl("Neg. ext.", sample_data(phyloseq$phyloseq_rel)$sample_type)
phyloseq_rel_nz <- subset_samples(phyloseq$phyloseq_rel, S.obs != 0)
#With all sampels
dec_f_all <- isContaminant(phyloseq_rel_nz, method="frequency", conc="DNA_bac_well")
dec_p_all <- isContaminant(phyloseq_rel_nz, method="prevalence", neg="is.neg", threshold=0.5)
dec_c_all <- isContaminant(phyloseq_rel_nz, method="combined", neg="is.neg", conc = "DNA_bac_well")
cat("decontam frequency - all sample")
## decontam frequency - all sample
dec_f_all %>% subset(.,.$contaminant)
cat("decontam prevalence - all sample")
## decontam prevalence - all sample
dec_p_all %>% subset(.,.$contaminant)
cat("decontam combined - all sample")
## decontam combined - all sample
dec_c_all %>% subset(.,.$contaminant)
#Stratified by sample type
cat("decontam prevalence - BAL")
## decontam prevalence - BAL
subset_samples(phyloseq_rel_nz, sample_type %in% c("BAL", "Neg. ext.")) %>%
isContaminant(., method="prevalence", neg = "is.neg", threshold = 0.5) %>% subset(.,.$contaminant)
cat("decontam prevalence - Nasal swab")
## decontam prevalence - Nasal swab
subset_samples(phyloseq_rel_nz, sample_type %in% c("Nasal swab", "Neg. ext.")) %>%
isContaminant(., method="prevalence", neg = "is.neg", threshold = 0.5) %>% subset(.,.$contaminant)
cat("decontam prevalence - Sputum")
## decontam prevalence - Sputum
subset_samples(phyloseq_rel_nz, sample_type %in% c("Sputum", "Neg. ext.")) %>%
isContaminant(., method="prevalence", neg = "is.neg", threshold = 0.5) %>% subset(.,.$contaminant)
cat("decontam frequency - BAL")
## decontam frequency - BAL
subset_samples(phyloseq_rel_nz, sample_type %in% c("BAL", "Neg. ext.")) %>%
isContaminant(method="frequency", conc="DNA_bac_well") %>% subset(.,.$contaminant)
cat("decontam frequency - Nasal swab")
## decontam frequency - Nasal swab
subset_samples(phyloseq_rel_nz, sample_type %in% c("Nasal swab", "Neg. ext.")) %>%
isContaminant(method="frequency", conc="DNA_bac_well") %>% subset(.,.$contaminant)
cat("decontam frequency - Sputum")
## decontam frequency - Sputum
subset_samples(phyloseq_rel_nz, sample_type %in% c("Sputum", "Neg. ext.")) %>%
isContaminant(method="frequency", conc="DNA_bac_well") %>% subset(.,.$contaminant)
cat("decontam combined - BAL")
## decontam combined - BAL
subset_samples(phyloseq_rel_nz, sample_type %in% c("BAL", "Neg. ext.")) %>%
isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant)
cat("decontam combined - Nasal swab")
## decontam combined - Nasal swab
subset_samples(phyloseq_rel_nz, sample_type %in% c("Nasal swab", "Neg. ext.")) %>%
isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant)
cat("decontam combined - Sputum")
## decontam combined - Sputum
subset_samples(phyloseq_rel_nz, sample_type %in% c("Sputum", "Neg. ext.")) %>%
isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant)
A7. LM of function alpha diversity
--> both stratified and nonstratified
sample_data <- sample_data(phyloseq$phyloseq_path_rpkm) %>% data.frame(check.names = F) %>% subset(., !is.nan(.$simpson))
cat("Species richness of all samples\n\n")
## Species richness of all samples
cat("S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample)\n\n")
## S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample)
lmer(S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|subject_id), data = sample_data) %>%
summary() %>%
.$coefficients %>%
data.frame(check.names = F) %>%
mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Estimate | Std. Error | df | t value | Pr(>|t|) | ||
|---|---|---|---|---|---|---|
| (Intercept) | -271.32533 | 63.94609 | 72.47916 | -4.2430322 | 0.0000643 |
|
| sample_typeNasal swab | 66.71179 | 39.44578 | 32.64336 | 1.6912273 | 0.1003209 | |
| sample_typeSputum | 142.80087 | 36.31835 | 41.80438 | 3.9319203 | 0.0003111 |
|
| treatmentlyPMA | 59.31768 | 31.18858 | 64.39962 | 1.9019036 | 0.0616583 | |
| treatmentBenzonase | 96.59262 | 31.65626 | 66.97650 | 3.0512965 | 0.0032643 |
|
| treatmentHost zero | 129.84121 | 32.25492 | 67.36194 | 4.0254694 | 0.0001466 |
|
| treatmentMolysis | 150.32258 | 32.65321 | 67.58822 | 4.6036081 | 0.0000189 |
|
| treatmentQIAamp | 86.27368 | 32.71224 | 67.61999 | 2.6373516 | 0.0103578 |
|
| log10(Final_reads) | 52.77944 | 11.01007 | 70.54314 | 4.7937427 | 0.0000088 |
|
| sample_typeNasal swab:treatmentlyPMA | -23.81209 | 39.44179 | 64.87664 | -0.6037273 | 0.5481301 | |
| sample_typeSputum:treatmentlyPMA | -5.02253 | 39.40167 | 62.83470 | -0.1274700 | 0.8989755 | |
| sample_typeNasal swab:treatmentBenzonase | -80.25929 | 38.01990 | 65.40988 | -2.1109811 | 0.0386006 |
|
| sample_typeSputum:treatmentBenzonase | -52.43718 | 38.65224 | 63.57018 | -1.3566400 | 0.1796940 | |
| sample_typeNasal swab:treatmentHost zero | -106.90271 | 36.81521 | 64.23327 | -2.9037647 | 0.0050492 |
|
| sample_typeSputum:treatmentHost zero | -76.35866 | 38.92443 | 62.75019 | -1.9617158 | 0.0542352 | |
| sample_typeNasal swab:treatmentMolysis | -106.97855 | 38.54815 | 66.00383 | -2.7751931 | 0.0071693 |
|
| sample_typeSputum:treatmentMolysis | -108.15918 | 39.34418 | 62.63915 | -2.7490513 | 0.0077999 |
|
| sample_typeNasal swab:treatmentQIAamp | -75.62038 | 36.68939 | 64.06591 | -2.0610968 | 0.0433576 |
|
| sample_typeSputum:treatmentQIAamp | -40.59820 | 38.59881 | 63.08680 | -1.0517993 | 0.2969039 |
#Species richness - stratified
cat("Species richness at NS\n\n")
## Species richness at NS
cat("S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(S.obs ~ treatment + log10 (Final_reads) + (1|subject_id), data = subset(sample_data, sample_data$sample_type == "Nasal swab")) %>%
summary() %>%
.$coefficients %>%
data.frame(check.names = F) %>%
mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Estimate | Std. Error | df | t value | Pr(>|t|) | ||
|---|---|---|---|---|---|---|
| (Intercept) | 68.35891 | 68.76371 | 28 | 0.9941131 | 0.3286817 | |
| treatmentlyPMA | 13.88909 | 17.11543 | 28 | 0.8114955 | 0.4239262 | |
| treatmentBenzonase | 20.95703 | 16.24183 | 28 | 1.2903126 | 0.2074965 | |
| treatmentHost zero | 56.43052 | 18.37320 | 28 | 3.0713495 | 0.0047048 |
|
| treatmentMolysis | 51.56131 | 16.32822 | 28 | 3.1578029 | 0.0037880 |
|
| treatmentQIAamp | 54.41631 | 19.71859 | 28 | 2.7596444 | 0.0100871 |
|
| log10(Final_reads) | 12.25116 | 10.45340 | 28 | 1.1719787 | 0.2510816 |
cat("Species richness at BAL\n\n")
## Species richness at BAL
cat("S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(S.obs ~ treatment + log10 (Final_reads) + (1|original_sample), data = subset(sample_data, sample_data$sample_type == "BAL")) %>%
summary() %>%
.$coefficients %>%
data.frame(check.names = F) %>%
mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Estimate | Std. Error | df | t value | Pr(>|t|) | ||
|---|---|---|---|---|---|---|
| (Intercept) | -640.83129 | 99.14027 | 19.00295 | -6.4638848 | 0.0000034 |
|
| treatmentlyPMA | 14.46191 | 33.76511 | 18.36963 | 0.4283092 | 0.6734060 | |
| treatmentBenzonase | 16.49562 | 35.65062 | 19.98741 | 0.4627021 | 0.6485757 | |
| treatmentHost zero | 39.81066 | 36.93641 | 19.98154 | 1.0778162 | 0.2939563 | |
| treatmentMolysis | 54.18626 | 37.78713 | 19.92262 | 1.4339874 | 0.1670859 | |
| treatmentQIAamp | -10.73947 | 37.91284 | 19.91125 | -0.2832674 | 0.7798961 | |
| log10(Final_reads) | 124.09731 | 17.94627 | 17.34430 | 6.9149360 | 0.0000022 |
|
cat("Species richness at sputum\n\n")
## Species richness at sputum
cat("S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## S.obs ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(S.obs ~ treatment + log10 (Final_reads) + (1|original_sample), data = subset(sample_data, sample_data$sample_type == "Sputum")) %>%
summary() %>%
.$coefficients %>%
data.frame(check.names = F) %>%
mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Estimate | Std. Error | df | t value | Pr(>|t|) | ||
|---|---|---|---|---|---|---|
| (Intercept) | -127.58303 | 198.98655 | 21.38568 | -0.6411641 | 0.5282282 | |
| treatmentlyPMA | 54.38241 | 28.77622 | 19.88708 | 1.8898385 | 0.0734460 | |
| treatmentBenzonase | 44.29211 | 36.29072 | 20.41549 | 1.2204804 | 0.2361930 | |
| treatmentHost zero | 53.75292 | 61.00733 | 20.95720 | 0.8810895 | 0.3882627 | |
| treatmentMolysis | 42.48496 | 71.13867 | 21.03201 | 0.5972132 | 0.5567405 | |
| treatmentQIAamp | 45.90423 | 52.95334 | 20.86328 | 0.8668807 | 0.3958683 | |
| log10(Final_reads) | 52.61786 | 33.96577 | 21.23460 | 1.5491436 | 0.1361234 |
#Shannon
cat("Shannon index\n\n")
## Shannon index
cat("Shannon ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## Shannon ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(data_shannon ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample), data = sample_data) %>%
summary() %>%
.$coefficients %>%
data.frame(check.names = F) %>%
mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Estimate | Std. Error | df | t value | Pr(>|t|) | ||
|---|---|---|---|---|---|---|
| (Intercept) | 0.6277442 | 0.1923582 | 72.90417 | 3.2634128 | 0.0016771 |
|
| sample_typeNasal swab | 0.3782276 | 0.0948640 | 66.24283 | 3.9870501 | 0.0001694 |
|
| sample_typeSputum | 0.3038993 | 0.0968430 | 65.80070 | 3.1380607 | 0.0025449 |
|
| treatmentlyPMA | 0.2080648 | 0.0912575 | 59.31888 | 2.2799751 | 0.0262196 |
|
| treatmentBenzonase | 0.1074536 | 0.0924809 | 65.10258 | 1.1619011 | 0.2495177 | |
| treatmentHost zero | 0.2663135 | 0.0943090 | 66.22861 | 2.8238390 | 0.0062621 |
|
| treatmentMolysis | 0.3493594 | 0.0955289 | 66.88889 | 3.6571063 | 0.0005034 |
|
| treatmentQIAamp | 0.0974529 | 0.0957099 | 66.98134 | 1.0182112 | 0.3122414 | |
| log10(Final_reads) | -0.0167122 | 0.0335618 | 71.90795 | -0.4979546 | 0.6200356 | |
| sample_typeNasal swab:treatmentlyPMA | -0.2199998 | 0.1167864 | 63.60638 | -1.8837795 | 0.0641674 | |
| sample_typeSputum:treatmentlyPMA | -0.0942524 | 0.1154524 | 55.89187 | -0.8163741 | 0.4177521 | |
| sample_typeNasal swab:treatmentBenzonase | -0.0695042 | 0.1120531 | 63.47768 | -0.6202792 | 0.5372946 | |
| sample_typeSputum:treatmentBenzonase | 0.0023288 | 0.1130661 | 57.07248 | 0.0205967 | 0.9836392 | |
| sample_typeNasal swab:treatmentHost zero | -0.2682060 | 0.1085902 | 60.90331 | -2.4698926 | 0.0163329 |
|
| sample_typeSputum:treatmentHost zero | -0.2089075 | 0.1142346 | 56.03885 | -1.8287583 | 0.0727592 | |
| sample_typeNasal swab:treatmentMolysis | -0.3849402 | 0.1134802 | 64.38350 | -3.3921345 | 0.0011901 |
|
| sample_typeSputum:treatmentMolysis | -0.2757981 | 0.1156639 | 56.27671 | -2.3844786 | 0.0204979 |
|
| sample_typeNasal swab:treatmentQIAamp | -0.1846564 | 0.1083181 | 60.79440 | -1.7047607 | 0.0933449 | |
| sample_typeSputum:treatmentQIAamp | -0.0636306 | 0.1130449 | 56.24157 | -0.5628795 | 0.5757540 |
#Simpson
cat("Inverse Simpson index\n\n")
## Inverse Simpson index
cat("InvSimp ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## InvSimp ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(data_invsimpson ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample), data = sample_data) %>%
summary() %>%
.$coefficients %>%
data.frame(check.names = F) %>%
mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Estimate | Std. Error | df | t value | Pr(>|t|) | ||
|---|---|---|---|---|---|---|
| (Intercept) | 2.5338586 | 0.3135385 | 70.96294 | 8.0814913 | 0.0000000 |
|
| sample_typeNasal swab | 0.2842938 | 0.1520494 | 68.03030 | 1.8697465 | 0.0658228 | |
| sample_typeSputum | 0.3782501 | 0.1558509 | 70.95435 | 2.4270003 | 0.0177657 |
|
| treatmentlyPMA | 0.3397324 | 0.1563898 | 61.65671 | 2.1723442 | 0.0336834 |
|
| treatmentBenzonase | 0.0472378 | 0.1563836 | 68.23317 | 0.3020637 | 0.7635223 | |
| treatmentHost zero | 0.2329645 | 0.1590151 | 69.37014 | 1.4650462 | 0.1474261 | |
| treatmentMolysis | 0.3567617 | 0.1607861 | 70.00663 | 2.2188598 | 0.0297400 |
|
| treatmentQIAamp | 0.0582910 | 0.1610497 | 70.09361 | 0.3619443 | 0.7184825 | |
| log10(Final_reads) | -0.1644187 | 0.0540196 | <<<<<<< HEAD 64.44675 ======= 64.44674 >>>>>>> a74da024aec8928028c246dbcfe35b96e9691229 | -3.0436866 | 0.0033807 |
|
| sample_typeNasal swab:treatmentlyPMA | -0.4423346 | 0.1983466 | 65.97618 | -2.2301099 | 0.0291475 |
|
| sample_typeSputum:treatmentlyPMA | -0.3759508 | 0.1993215 | 57.44567 | -1.8861530 | 0.0643330 | |
| sample_typeNasal swab:treatmentBenzonase | -0.0381545 | 0.1903215 | 66.04230 | -0.2004738 | 0.8417261 | |
| sample_typeSputum:treatmentBenzonase | -0.1934003 | 0.1947030 | 58.93845 | -0.9933096 | 0.3246198 | |
| sample_typeNasal swab:treatmentHost zero | -0.3256749 | 0.1855479 | 62.96653 | -1.7552062 | 0.0840880 | |
| sample_typeSputum:treatmentHost zero | -0.4518719 | 0.1971485 | 57.67798 | -2.2920381 | 0.0255722 |
|
| sample_typeNasal swab:treatmentMolysis | -0.5649464 | 0.1923302 | 67.03548 | -2.9373781 | 0.0045336 |
|
| sample_typeSputum:treatmentMolysis | -0.5317988 | 0.1995233 | 57.91638 | -2.6653467 | 0.0099491 |
|
| sample_typeNasal swab:treatmentQIAamp | -0.1812611 | 0.1851750 | 62.57751 | -0.9788638 | 0.3314165 | |
| sample_typeSputum:treatmentQIAamp | -0.3137393 | 0.1950063 | 57.95487 | -1.6088678 | 0.1130797 |
#BPI
cat("Beger Parker index\n\n")
## Beger Parker index
cat("BPI ~ sample_type + log10 (Final_reads) + (1|original_sample)\n\n")
## BPI ~ sample_type + log10 (Final_reads) + (1|original_sample)
lmer(dbp ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample), data = sample_data) %>%
summary() %>%
.$coefficients %>%
data.frame(check.names = F) %>%
mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Estimate | Std. Error | df | t value | Pr(>|t|) | ||
|---|---|---|---|---|---|---|
| (Intercept) | 0.3440391 | 0.1258762 | 70.13993 | 2.7331535 | 0.0079332 |
|
| sample_typeNasal swab | -0.0940641 | 0.0609175 | 68.34780 | -1.5441228 | 0.1271764 | |
| sample_typeSputum | -0.1495148 | 0.0626267 | 71.73071 | -2.3873954 | 0.0196055 |
|
| treatmentlyPMA | -0.1360672 | 0.0637784 | 62.70891 | -2.1334365 | 0.0368082 |
|
| treatmentBenzonase | -0.0137189 | 0.0634994 | 69.15451 | -0.2160472 | 0.8295870 | |
| treatmentHost zero | -0.0775264 | 0.0645073 | 70.21212 | -1.2018236 | 0.2334698 | |
| treatmentMolysis | -0.1514704 | 0.0651879 | 70.79134 | -2.3235979 | 0.0230235 |
|
| treatmentQIAamp | -0.0147607 | 0.0652893 | 70.86956 | -0.2260821 | 0.8217883 | |
| log10(Final_reads) | 0.0701153 | 0.0215697 | 61.84125 | 3.2506433 | 0.0018652 |
|
| sample_typeNasal swab:treatmentlyPMA | 0.1871812 | 0.0806756 | 66.80034 | 2.3201722 | 0.0233951 |
|
| sample_typeSputum:treatmentlyPMA | 0.1784389 | 0.0814836 | 58.41893 | 2.1898742 | 0.0325394 |
|
| sample_typeNasal swab:treatmentBenzonase | 0.0181788 | 0.0774064 | 66.92125 | 0.2348489 | 0.8150431 | |
| sample_typeSputum:treatmentBenzonase | 0.0966450 | 0.0795283 | 59.95203 | 1.2152280 | 0.2290433 | |
| sample_typeNasal swab:treatmentHost zero | 0.1112523 | 0.0756141 | 63.85620 | 1.4713183 | 0.1461162 | |
| sample_typeSputum:treatmentHost zero | 0.1766839 | 0.0805844 | 58.67215 | 2.1925319 | 0.0323198 |
|
| sample_typeNasal swab:treatmentMolysis | 0.2327702 | 0.0781693 | 67.88223 | 2.9777684 | 0.0040230 |
|
| sample_typeSputum:treatmentMolysis | 0.2309710 | 0.0815446 | 58.89704 | 2.8324495 | 0.0063138 |
|
| sample_typeNasal swab:treatmentQIAamp | 0.0548290 | 0.0754814 | 63.39830 | 0.7263905 | 0.4702729 | |
| sample_typeSputum:treatmentQIAamp | 0.1291496 | 0.0796962 | 58.96349 | 1.6205246 | 0.1104560 |
A8. permanova of function beta diversity
--> both stratified and nonstratified
phyloseq_rel_nz <- subset_samples(phyloseq$phyloseq_path_rpkm, S.obs != 0 & sample_type %in% c("BAL", "Nasal swab", "Sputum"))
bray_perm_uni <- vegan::adonis2(distance(phyloseq_rel_nz, method="bray") ~ sample_type + log10(Final_reads) + treatment + subject_id,
data = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F), permutations = 10000)
bray_perm <- vegan::adonis2(distance(phyloseq_rel_nz, method="bray") ~ sample_type + log10(Final_reads) + lypma + benzonase + host_zero + molysis + qiaamp + subject_id,
data = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F), permutations = 10000)
bray_perm_strata <- vegan::adonis2(distance(phyloseq_rel_nz, method="bray") ~ sample_type + log10(Final_reads) + lypma + benzonase + host_zero + molysis + qiaamp,
data = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F),
strata = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F) %>% .$subject_id, permutations = 10000)
bray_perm_inter <- vegan::adonis2(distance(phyloseq_rel_nz, method="bray") ~ sample_type * treatment + log10(Final_reads) + subject_id,
data = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F), permutations = 10000)
bray_perm_ns <- vegan::adonis2(distance(subset_samples(phyloseq_rel_nz, sample_type == "Nasal swab"), method="bray") ~ lypma + benzonase + host_zero + molysis + qiaamp + log10(Final_reads) + subject_id,
data = subset_samples(phyloseq_rel_nz, sample_type == "Nasal swab") %>% sample_data %>% data.frame(check.names = F), permutations = 10000)
bray_perm_bal <- vegan::adonis2(distance(subset_samples(phyloseq_rel_nz, sample_type == "BAL"), method="bray") ~ lypma + benzonase + host_zero + molysis + qiaamp + log10(Final_reads) + subject_id,
data = subset_samples(phyloseq_rel_nz, sample_type == "BAL") %>% sample_data %>% data.frame(check.names = F), permutations = 10000)
bray_perm_spt <- vegan::adonis2(distance(subset_samples(phyloseq_rel_nz, sample_type == "Sputum"), method="bray") ~ lypma + benzonase + host_zero + molysis + qiaamp + log10(Final_reads) + subject_id,
data = subset_samples(phyloseq_rel_nz, sample_type == "Sputum") %>% sample_data %>% data.frame(check.names = F), permutations = 10000)
cat("\nUnivariate analysis\n")
##
## Univariate analysis
bray_perm_uni %>% data.frame(check.names = F) %>% rownames_to_column("row.names") %>%
mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
row.names == "treatment" ~ 'Treatment',
row.names == "subject_id" ~ 'Subject',
row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
row.names == "Residual" ~ 'Residual',
row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>%
round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Df | SumOfSqs | R2 | F | Pr(>F) | ||
|---|---|---|---|---|---|---|
| Sample type | 2 | 1.557 | 0.344 | 53.603 | 0.00 |
|
| log10(Final reads) | 1 | 1.024 | 0.226 | 70.491 | 0.00 |
|
| Treatment | 5 | 0.124 | 0.027 | 1.703 | 0.14 | |
| Subject | 10 | 0.766 | 0.169 | 5.276 | 0.00 |
|
| Residual | 73 | 1.060 | 0.234 | NA | NA | |
| Total | 91 | 4.530 | 1.000 | NA | NA |
cat("\nDetailed treatment\n")
##
## Detailed treatment
bray_perm %>% data.frame(check.names = F) %>% rownames_to_column("row.names") %>%
mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
row.names == "lypma" ~ 'lyPMA',
row.names == "benzonase" ~ 'Benzonase',
row.names == "host_zero" ~ 'Host zero',
row.names == "molysis" ~ 'Molysis',
row.names == "qiaamp" ~ 'QIAamp',
row.names == "subject_id" ~ 'Subject',
row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
row.names == "Residual" ~ 'Residual',
row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>%
round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Df | SumOfSqs | R2 | F | Pr(>F) | ||
|---|---|---|---|---|---|---|
| Sample type | 2 | 1.557 | 0.344 | 53.603 | 0.000 |
|
| log10(Final reads) | 1 | 1.024 | 0.226 | 70.491 | 0.000 |
|
| lyPMA | 1 | 0.019 | 0.004 | 1.284 | 0.263 | |
| Benzonase | 1 | 0.008 | 0.002 | 0.579 | 0.459 | |
| Host zero | 1 | 0.008 | 0.002 | 0.564 | 0.460 | |
| Molysis | 1 | 0.087 | 0.019 | 6.020 | 0.015 |
|
| QIAamp | 1 | 0.001 | 0.000 | 0.067 | 0.836 | |
| Subject | 10 | 0.766 | 0.169 | 5.276 | 0.000 |
|
| Residual | 73 | 1.060 | 0.234 | NA | NA | |
| Total | 91 | 4.530 | 1.000 | NA | NA |
cat("\n Strata -detailed treatment\n")
##
## Strata -detailed treatment
bray_perm_strata %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>%
mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
row.names == "lypma" ~ 'lyPMA',
row.names == "benzonase" ~ 'Benzonase',
row.names == "host_zero" ~ 'Host zero',
row.names == "molysis" ~ 'Molysis',
row.names == "qiaamp" ~ 'QIAamp',
row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
row.names == "Residual" ~ 'Residual',
row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>%
round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Df | SumOfSqs | R2 | F | Pr(>F) | ||
|---|---|---|---|---|---|---|
| Sample type | 2 | 1.557 | 0.344 | 35.378 | 0.000 |
|
| log10(Final reads) | 1 | 1.024 | 0.226 | 46.525 | 0.000 |
|
| lyPMA | 1 | 0.019 | 0.004 | 0.847 | 0.424 | |
| Benzonase | 1 | 0.008 | 0.002 | 0.382 | 0.521 | |
| Host zero | 1 | 0.008 | 0.002 | 0.372 | 0.570 | |
| Molysis | 1 | 0.087 | 0.019 | 3.973 | 0.054 | |
| QIAamp | 1 | 0.001 | 0.000 | 0.044 | 0.964 | |
| Residual | 83 | 1.826 | 0.403 | NA | NA | |
| Total | 91 | 4.530 | 1.000 | NA | NA |
cat("\n Interaction term \n")
##
## Interaction term
bray_perm_inter %>% data.frame(check.names = F) %>% rownames_to_column("row.names") %>%
mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
row.names == "treatment" ~ 'Treatment',
row.names == "subject_id" ~ 'Subject',
row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
row.names == "sample_type:treatment" ~ 'Sample type X treatment',
row.names == "Residual" ~ 'Residual',
row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>%
round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Df | SumOfSqs | R2 | F | Pr(>F) | ||
|---|---|---|---|---|---|---|
| Sample type | 2 | 1.557 | 0.344 | 66.198 | 0.000 |
|
| Treatment | 5 | 0.684 | 0.151 | 11.631 | 0.000 |
|
| log10(Final reads) | 1 | 0.463 | 0.102 | 39.412 | 0.000 |
|
| Subject | 10 | 0.766 | 0.169 | 6.515 | 0.000 |
|
| Sample type X treatment | 10 | 0.319 | 0.070 | 2.715 | 0.008 |
|
| Residual | 63 | 0.741 | 0.164 | NA | NA | |
| Total | 91 | 4.530 | 1.000 | NA | NA |
cat("\n Stratified - nasal swab \n")
##
## Stratified - nasal swab
bray_perm_ns %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>%
mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
row.names == "benzonase" ~ 'Benzonase',
row.names == "host_zero" ~ 'Host zero',
row.names == "molysis" ~ 'Molysis',
row.names == "qiaamp" ~ 'QIAamp',
row.names == "subject_id" ~ 'Subject id',
row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
row.names == "Residual" ~ 'Residual',
row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>%
round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Df | SumOfSqs | R2 | F | Pr(>F) | ||
|---|---|---|---|---|---|---|
| lyPMA | 1 | 0.011 | 0.042 | 2.049 | 0.147 | |
| Benzonase | 1 | 0.011 | 0.040 | 1.940 | 0.159 | |
| Host zero | 1 | 0.009 | 0.033 | 1.615 | 0.205 | |
| Molysis | 1 | 0.014 | 0.052 | 2.566 | 0.104 | |
| QIAamp | 1 | 0.054 | 0.197 | 9.648 | 0.003 |
|
| log10(Final reads) | 1 | 0.028 | 0.103 | 5.051 | 0.027 |
|
| Subject id | 2 | 0.001 | 0.003 | 0.083 | 0.980 | |
| Residual | 26 | 0.146 | 0.530 | NA | NA | |
| Total | 34 | 0.275 | 1.000 | NA | NA |
cat("\n Stratified - BAL \n")
##
## Stratified - BAL
bray_perm_bal %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>%
mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
row.names == "benzonase" ~ 'Benzonase',
row.names == "host_zero" ~ 'Host zero',
row.names == "molysis" ~ 'Molysis',
row.names == "qiaamp" ~ 'QIAamp',
row.names == "subject_id" ~ 'Subject id',
row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
row.names == "Residual" ~ 'Residual',
row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>%
round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Df | SumOfSqs | R2 | F | Pr(>F) | ||
|---|---|---|---|---|---|---|
| lyPMA | 1 | 0.040 | 0.020 | 1.305 | 0.268 | |
| Benzonase | 1 | 0.115 | 0.058 | 3.716 | 0.070 | |
| Host zero | 1 | 0.027 | 0.014 | 0.886 | 0.370 | |
| Molysis | 1 | 0.180 | 0.091 | 5.820 | 0.028 |
|
| QIAamp | 1 | 0.001 | 0.000 | 0.029 | 0.904 | |
| log10(Final reads) | 1 | 0.603 | 0.306 | 19.535 | 0.001 |
|
| Subject id | 4 | 0.513 | 0.260 | 4.153 | 0.017 |
|
| Residual | 16 | 0.494 | 0.250 | NA | NA | |
| Total | 26 | 1.973 | 1.000 | NA | NA |
cat("\n Stratified - sputum \n")
##
## Stratified - sputum
bray_perm_spt %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>%
mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
row.names == "benzonase" ~ 'Benzonase',
row.names == "host_zero" ~ 'Host zero',
row.names == "molysis" ~ 'Molysis',
row.names == "qiaamp" ~ 'QIAamp',
row.names == "subject_id" ~ 'Subject id',
row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
row.names == "Residual" ~ 'Residual',
row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>%
round(3) %>% mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
.default = " ")) %>%
kbl(format = "html") %>%
kable_styling(full_width = 0, html_font = "serif")
| Df | SumOfSqs | R2 | F | Pr(>F) | ||
|---|---|---|---|---|---|---|
| lyPMA | 1 | 0.017 | 0.023 | 3.920 | 0.053 | |
| Benzonase | 1 | 0.002 | 0.003 | 0.440 | 0.557 | |
| Host zero | 1 | 0.064 | 0.089 | 15.072 | 0.001 |
|
| Molysis | 1 | 0.139 | 0.191 | 32.510 | 0.000 |
|
| QIAamp | 1 | 0.370 | 0.509 | 86.595 | 0.000 |
|
| log10(Final reads) | 1 | 0.032 | 0.045 | 7.588 | 0.011 |
|
| Subject id | 4 | 0.021 | 0.029 | 1.224 | 0.332 | |
| Residual | 19 | 0.081 | 0.112 | NA | NA | |
| Total | 29 | 0.726 | 1.000 | NA | NA |
A9. DA analysis for function, by sample type and treatment
--> both stratified and nonstratified
--> multiple levels?
No normalization will be conducted for RPKM data.
Thread: https://forum.biobakery.org/t/maaslin-with-shortbred-results-and-panphlan/3102/2
<<<<<<< HEADDone.
#===============================================================================
#BTC.LineZero.Footer.1.1.0
#===============================================================================
#R markdown citation generator.
#===============================================================================
#RLB.Dependencies:
# magrittr, pacman, stringr
#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
#BTC.Dependencies:
# LineZero.Header
#===============================================================================
#Generates citations for each explicitly loaded library.
#=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
str_libraries <- c("r", str_libraries)
for (str_libraries in str_libraries) {
str_libraries |>
pacman::p_citation() |>
print(bibtex = FALSE) |>
capture.output() %>%
.[-1:-3] %>% .[. != ""] |>
stringr::str_squish() |>
stringr::str_replace("_", "") |>
cat()
cat("\n")
}
## R Core Team (2022). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/. We have invested a lot of time and effort in creating R, please cite it when using it for data analysis. See also 'citation("pkgname")' for citing R packages.
## Wickham H, Bryan J (2023). readxl: Read Excel Files_. R package version 1.4.2, <https://CRAN.R-project.org/package=readxl>.
## phyloseq: An R package for reproducible interactive analysis and graphics of microbiome census data. Paul J. McMurdie and Susan Holmes (2013) PLoS ONE 8(4):e61217.
## Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). "Welcome to the tidyverse." Journal of Open Source Software_, *4*(43), 1686. doi:10.21105/joss.01686 <https://doi.org/10.21105/joss.01686>.
## Rinker, T. W. & Kurkiewicz, D. (2017). pacman: Package Management for R. version 0.5.0. Buffalo, New York. http://github.com/trinker/pacman
## Garbett SP, Stephens J, Simonov K, Xie Y, Dong Z, Wickham H, Horner J, reikoch, Beasley W, O'Connor B, Warnes GR, Quinn M, Kamvar ZN (2023). yaml: Methods to Convert R Data to YAML and Back_. R package version 2.3.7, <https://CRAN.R-project.org/package=yaml>. ATTENTION: This citation information has been auto-generated from the package DESCRIPTION file and may need manual editing, see 'help("citation")'.
## H. Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016.
## Oksanen J, Simpson G, Blanchet F, Kindt R, Legendre P, Minchin P, O'Hara R, Solymos P, Stevens M, Szoecs E, Wagner H, Barbour M, Bedward M, Bolker B, Borcard D, Carvalho G, Chirico M, De Caceres M, Durand S, Evangelista H, FitzJohn R, Friendly M, Furneaux B, Hannigan G, Hill M, Lahti L, McGlinn D, Ouellette M, Ribeiro Cunha E, Smith T, Stier A, Ter Braak C, Weedon J (2022). vegan: Community Ecology Package. R package version 2.6-4, <https://CRAN.R-project.org/package=vegan>.
## Leo Lahti et al. microbiome R package. URL: http://microbiome.github.io
## Kassambara A (2023). ggpubr: 'ggplot2' Based Publication Ready Plots. R package version 0.6.0, <https://CRAN.R-project.org/package=ggpubr>.
## Simon Garnier, Noam Ross, Robert Rudis, Antônio P. Camargo, Marco Sciaini, and Cédric Scherer (2021). Rvision - Colorblind-Friendly Color Maps for R. R package version 0.6.2.
## Davis NM, Proctor D, Holmes SP, Relman DA, Callahan BJ (2017). "Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data." bioRxiv_, 221499. doi:10.1101/221499 <https://doi.org/10.1101/221499>.
## Auguie B (2017). gridExtra: Miscellaneous Functions for "Grid" Graphics. R package version 2.3, <https://CRAN.R-project.org/package=gridExtra>.
## Kassambara A (2023). ggpubr: 'ggplot2' Based Publication Ready Plots. R package version 0.6.0, <https://CRAN.R-project.org/package=ggpubr>.
## Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1-48. doi:10.18637/jss.v067.i01.
## Kuznetsova A, Brockhoff PB, Christensen RHB (2017). "lmerTest Package: Tests in Linear Mixed Effects Models." Journal of Statistical Software, *82*(13), 1-26. doi:10.18637/jss.v082.i13 <https://doi.org/10.18637/jss.v082.i13>.
## Ooms J (2023). writexl: Export Data Frames to Excel 'xlsx' Format_. R package version 1.4.2, <https://CRAN.R-project.org/package=writexl>.
## Gonçalves da Silva A (2017). harrietr: Wrangle Phylogenetic Distance Matrices and Other Utilities. R package version 0.2.3, <https://CRAN.R-project.org/package=harrietr>.
## Mallick H et al. (2020). Multivariable Association in Population-scale Meta-omics Studies, http://huttenhower.sph.harvard.edu/maaslin2. To cite the MaAsLin 2 software, please use: Mallick H, Rahnavard A, McIver LJ (2020). MaAsLin 2: Multivariable Association in Population-scale Meta-omics Studies. R/Bioconductor package, http://huttenhower.sph.harvard.edu/maaslin2.
## Wilke C, Wiernik B (2022). ggtext: Improved Text Rendering Support for 'ggplot2'. R package version 0.1.2, <https://CRAN.R-project.org/package=ggtext>.
## Aphalo P (2022). ggpmisc: Miscellaneous Extensions to 'ggplot2'_. R package version 0.5.2, <https://CRAN.R-project.org/package=ggpmisc>.
## Auguie B (2017). gridExtra: Miscellaneous Functions for "Grid" Graphics. R package version 2.3, <https://CRAN.R-project.org/package=gridExtra>.
## Wood S, Scheipl F (2020). gamm4: Generalized Additive Mixed Models using 'mgcv' and 'lme4'. R package version 0.2-6, <https://CRAN.R-project.org/package=gamm4>. ATTENTION: This citation information has been auto-generated from the package DESCRIPTION file and may need manual editing, see 'help("citation")'.
## Hadley Wickham (2007). Reshaping Data with the reshape Package. Journal of Statistical Software, 21(12), 1-20. URL http://www.jstatsoft.org/v21/i12/.
## Zhu H (2021). kableExtra: Construct Complex Table with 'kable' and Pipe Syntax. R package version 1.3.4, <https://CRAN.R-project.org/package=kableExtra>.
## Yihui Xie (2023). knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.42. Yihui Xie (2015) Dynamic Documents with R and knitr. 2nd edition. Chapman and Hall/CRC. ISBN 978-1498716963 Yihui Xie (2014) knitr: A Comprehensive Tool for Reproducible Research in R. In Victoria Stodden, Friedrich Leisch and Roger D. Peng, editors, Implementing Reproducible Computational Research. Chapman and Hall/CRC. ISBN 978-1466561595
#===============================================================================